-
Notifications
You must be signed in to change notification settings - Fork 0
/
testspec.f90
232 lines (232 loc) · 9.82 KB
/
testspec.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
module specfiles
integer(kind=4) :: iunit,iscan,ispecerr,ncolumns
integer(kind=4),parameter :: LINELENGTH=512, WORDLENGTH=100
integer(kind=4),parameter :: NWORDS=60, NLINES=50
integer(kind=4),dimension(NLINES) :: headerwords
real(kind=8) :: scanstart, scanend
real(kind=8),dimension(NWORDS,NLINES) :: wordvalues
real(kind=8),dimension(NWORDS) :: Q
character(len=6) :: lnfmt
character(len=LINELENGTH) :: filnam, line
character(len=LINELENGTH) :: specsfilename, specsdate
character(len=WORDLENGTH) :: scantype
character(len=WORDLENGTH),dimension(NWORDS,NLINES) :: words
character(len=WORDLENGTH),dimension(NWORDS) :: columnlabels
integer(kind=4) :: i0, i9, id, isp, im
character(len=WORDLENGTH) :: FIRSTDET,LASTDET,TWOTTH
character(len=WORDLENGTH) :: MONITORCOL="Monitor"
logical topofscan
data lnfmt /'(a512)'/
data iunit /15/
contains
subroutine getfile
logical :: od
! See if there is already a file open, if there is then close it
inquire(unit=iunit,opened=od)
if(od)close(iunit)
! filnam must already be filled in, or we'll try to open garbage
1 open(unit=iunit,file=filnam,form='formatted', &
& access='sequential', status='old',iostat=ispecerr)
if(ispecerr.ne.0) then
write(*,'(3a,i6)')'Problem with file ',filnam(1:len_trim(filnam)),&
& ' iostat=',ispecerr
write(*,'(a)')'Please supply a valid SPEC file name'
read(*,lnfmt)filnam
goto 1
endif
topofscan=.false.
! Supply defaults in case dumb specfile has no #L
ncolumns=14 ; columnlabels(1)='2_theta'; columnlabels(2)='Epoch'
columnlabels(3)='Seconds'; columnlabels(4)='MA0'
columnlabels(5)='MA1'; columnlabels(6)='MA2'
columnlabels(7)='MA3'; columnlabels(8)='MA4'
columnlabels(9)='MA5'; columnlabels(10)='MA6'
columnlabels(11)='MA7'; columnlabels(12)='MA8'
columnlabels(13)='Monitor'; columnlabels(14)='Fluo det'
i0=ichar('0');i9=ichar('9')
id=ichar('.');isp=ichar(' ');im=ichar('-')
return
end subroutine getfile
subroutine findscan(n)
integer(kind=4),intent(in)::n
character(len=7) :: rd ! enough space for yes, no and unknown
inquire(unit=iunit,read=rd) ! Can we read the file?
if(rd(1:3).ne.'YES') call getfile ! If not go open it
if(iscan.eq.n .and. topofscan) return
if(iscan.gt.n) rewind(n) !
2 call readheader ! in case top of file
if(iscan.ne.n)then ! work through file to find #
1 read(iunit,lnfmt,err=10,end=10)line
if(line(1:1).eq.'#') goto 2
go to 1
else
! Scan found - check scantype elsewhere - not a job for specfile module
return
endif
10 write(*,*)
write(*,'(a,i5,a)')'Scan ',n,' not found in file '// &
& filnam(1:len_trim(filnam))
ispecerr=-2
return
end subroutine findscan
subroutine getdata(a,n)
! A is an array, dimensioned by the number of data items to be expected
integer(kind=4),intent(in) :: n
real(kind=8),dimension(n),intent(out) :: a
integer(kind=4), parameter :: one = 1
if(topofscan)then; topofscan=.false. ; else
read(iunit,lnfmt,err=10,end=10)line ! goes via line
topofscan=.false.
endif
call rdnums(one,n,a) ! reads from line
if(line(1:1).eq.'#')then; ispecerr=-1 ;goto 100 ; endif
go to 100
10 ispecerr=-1 ! End of file
100 return ; end subroutine getdata
subroutine readheader
character(len=1) :: letter
character(len=128) :: motor
integer(kind=4) :: i,j
real(kind=8) :: a
integer(kind=4),parameter :: four=4
if(line(1:1).eq.'#') goto 2 ! if already on header, don't skip
1 read(iunit,lnfmt,end=100)line
2 if(line(1:1).eq.'#') then; letter=line(2:2); topofscan=.true.
select case((letter))
case('O')
read(line(3:3),'(i1)')i
call split(line(4:LINELENGTH),words(:,i+1),WORDLENGTH,NWORDS,j)
headerwords(i+1)=j
case('P')
read(line(3:3),'(i1)')i
call rdnums(four,NWORDS,wordvalues(:,i+1))
case('Q')
! Copy the Q from the header into our specfile module
read(line(3:LINELENGTH),*,end=1, err=1) Q
case('N')
read(line(3:len_trim(line)),*,end=1)ncolumns
case('L') ! signals end of header, bug out here !
call split(line(3:LINELENGTH),columnlabels,WORDLENGTH,NWORDS,i)
if(ncolumns.ne.i) &
& write(*,'(a,i5)')'error reading header for scan',iscan
case('S')
read(line(3:len_trim(line)),*,end=1)iscan,scantype,motor, &
& scanstart, scanend
case('F')
specsfilename=line(3:len_trim(line)) ! get original filename
case('D')
read(line(3:len_trim(line)),'(a256)',end=1)specsdate
case('C'); continue ! comments
case('G'); continue ! no idea - always zero
case('E'); continue ! epoch - we don't care what it was for now.
case default
end select
else ! Not a # line, so assume end of header
read(line,*,err=1,end=1)a ! catch blank lines in header
topofscan=.true. !!! Must be able to read a number !
return !!! escapes from routine here !!!!!
endif
goto 1
100 ispecerr=-1; return; end subroutine readheader
real(kind=8) function getheadervalue(string)
character(len=*),intent(in) :: string
integer(kind=4) :: i, j, k
k=len(string); getheadervalue=0.0
do i=1,NWORDS; do j=1,NLINES
if(string(1:k).eq.words(i,j)(1:k))then
getheadervalue=wordvalues(i,j); return
endif
enddo; enddo
return; end function getheadervalue
integer(kind=4) function whichcolumn(string)
! Interprets the #L line information
character(len=*),intent(in) :: string
integer(kind=4) :: i, j
j=len(string) ; whichcolumn=-1
do i=1,NWORDS
if(string(1:j).eq.columnlabels(i)(1:j))then
whichcolumn=i ; return
endif
enddo
return; end function whichcolumn
subroutine rdnums(ic,n,values)
! Placed here in a subroutine in case of formatting or error handling problems
integer(kind=4),intent(in) :: n, ic
real(kind=8),dimension(n),intent(inout) :: values
if(line(ic:ic).eq.'#')then; ispecerr=-1; return; endif
read(line(ic:len_trim(line)),*,err=10,end=20)values(1:ncolumns)
goto 100
10 write(*,*)'Error reading line, looking for ',ncolumns,' values'
write(*,*)line(ic:len_trim(line))
write(*,*)values
20 continue
100 return
end subroutine rdnums
subroutine split(instring,outstrings,lenout,n,i)
integer(kind=4),intent(in) :: lenout,n
character(len=*),intent(in) :: instring
character(len=lenout),dimension(n),intent(out) :: outstrings
integer(kind=4),intent(out) :: i
integer(kind=4) :: j,k,l
j=1; k=1
do i=1,len_trim(instring) ! hope len > len_trim or array overstep
if(instring(i:i+1).ne.' ')then
outstrings(j)(k:k)=instring(i:i)
if(k.ne.1 .and. k.lt.lenout) k=k+1
if(instring(i:i).ne.' '.and.k.eq.1)k=k+1
else
if(k.gt.1)then ; do l=k,lenout
outstrings(j)(l:l)=' ' ! blank pad end of string
enddo ; j=j+1; k=1; if(j.gt.n)exit ; endif
endif
enddo
! Blank pad last string if necessary
if(k.gt.1)then
do l=k,lenout;outstrings(j)(l:l)=' ';enddo
endif
i=j; return; end subroutine split
end module specfiles
program testspec
use specfiles
integer(kind=4) :: i,j
real(kind=8),allocatable :: data(:), prevdata(:)
real :: time1, time2
call getarg(1,filnam)
call getfile
call cpu_time(time1)
call readheader ! Get's the header at the top of the file
2 write(*,'(a5,i5)')'Scan ',iscan
write(*,'(a)')'Header words and values'
do j=1,NLINES
if(headerwords(j).gt.0)then
do i=1,headerwords(j)
write(*,'(a20,f16.8)')words(i,j),wordvalues(i,j)
enddo
endif ! headerwords(j.gt.0)
enddo
write(*,'(a,i5)')'Number of columns = ', ncolumns
if(.not.allocated(data))allocate(data(ncolumns))
if(.not.allocated(prevdata))allocate(prevdata(ncolumns))
if(size(data).ne.ncolumns)then
deallocate(data); deallocate(prevdata)
allocate(data(ncolumns));allocate(prevdata(ncolumns))
endif
j=0
1 prevdata=data ! Start of loop through reading data
call getdata(data,ncolumns)
j=j+1
if(ispecerr.eq.0)goto 1 ! loop to end of scan
ispecerr=0 ! OK, we reached the end of the scan
write(*,'(a)')'Last data line in scan'
do i=1,ncolumns
write(*,'(a20,f16.8)')columnlabels(i),prevdata(i)
enddo
write(*,'(a,i10)')'Number of data points = ',j
call findscan(iscan+1) ! get the next scan (assumes sequential)
if(ispecerr.eq.0)goto 2 ! process next scan
deallocate(data); deallocate(prevdata)
call cpu_time(time2)
write(*,*)
write(*,'(a,f8.4,a)')'Time to read file = ',time2-time1,' seconds'
write(*,'(a)')'Thats all folks!'
end program testspec