Files
fit/gen/dat_sics.f
2022-08-19 15:22:33 +02:00

318 lines
7.1 KiB
Fortran

subroutine dat_sics
c -------------------
external dat_sics_desc
external dat_sics_opts
external dat_get_datanumber
external dat_sics_read
integer dtype/0/
call dat_init_desc(dtype, dat_sics_desc)
call dat_init_opts(dtype, dat_sics_opts)
call dat_init_high(dtype, dat_get_datanumber)
call dat_init_read(dtype, dat_sics_read)
end
subroutine dat_sics_desc(text)
! ------------------------------
character*(*) text ! (out) description
! type description
! ----------------------------------
text='SICS SICS-ASCII (TOPSI,TriCS)'
end
subroutine dat_sics_opts
! ------------------------
print '(x,a)'
1,'x,y: x-axis and y-axis column name'
end
subroutine dat_sics_read
1 (lun, forced, nread, putval, nmax, xx, yy, ss, ww)
! ----------------------------------------------------
implicit none
integer lun ! (in) logical unit number (file will be closed if successful)
integer forced ! 0: read only if type is sure; 1: forced read
integer nread ! (out) >=0: = number of points read, file closed
! -1: not correct type, file rewinded
! -2: correct type, but unreadable, file rewinded
external putval ! (in) subroutine to put name/value pairs.
! for numeric data: call putval('name', value) ! value must be real
! for character data: call putval('name=text', 0.0)
integer nmax ! max. number of points
real xx(*) ! x-values
real yy(*) ! y-values
real ss(*) ! sigma
real ww(*) ! weights (original monitor)
! local
integer mcol
parameter (mcol=32)
real y,s,r,f,ymon,values(mcol),val
integer i,j,l,errcnt,ncol,ycol,pcol,xcol,ycol2,ipol,xcol0
integer iostat, lz
logical pol
character line*132, preset*16, col2*16
character xaxis*16, xaxis0*16, yaxis*16, yaxis2*16
read(lun,'(a)',err=100,end=100) line
if (line(1:16) .eq. '##SICS ASCII at ') then
call putval('Instrument='//line(17:),0.0)
else
i=index(line,'Data File ****')
if (i .eq. 0) i=index(line,'SCAN File ****')
if (i .eq. 0) then
if (forced .le. 0) goto 100
else
j=index(line,'*** ')
if (j .lt. i) call putval('Instrument='//line(j+4:i-2),0.0)
endif
endif
nread=0
errcnt=0
xcol=0
xcol0=0
1 read(lun, '(a)', err=99,end=99) line
iostat=1
if (line(1:20) .ne. 'Scanning Variables: '
1 .and. line(1:20) .ne. 'scanning variables: ') then
i=index(line,'=')
if (i .le. 1) goto 1
call str_first_nonblank(line(i+1:), j)
call str_trim(line(1:i-1), line(1:i-1), l)
iostat=1
if (line(1:l) .eq. 'Sample Name') then
l=6
elseif (line(1:l) .eq. 'Original Filename' .or.
1 line(1:l) .eq. 'original_filename') then
goto 1
elseif (line(1:l) .eq. 'Title' .or.
1 line(1:l) .eq. 'title') then
call dat_group(1, putval)
elseif (line(1:13) .eq. 'File Creation' .or.
1 line(1:4) .eq. 'date') then
line(1:l)='Date'
l=4
else if (line(i+j:i+j) .eq. '-' .or.
1 line(i+j:i+j) .ge. '0' .and.
1 line(i+j:i+j) .le. '9') then
if (line(1:l) .eq. 'Sample Theta') then
line(1:l)='2-theta'
l=7
else if (line(1:l) .eq. 'Temperature' .or.
1 line(1:l) .eq. 'temp') then
line(1:l)='Temp'
l=4
endif
lz=index(line(1:l),' ')
if (lz .ne. 0) then
if (line(lz:lz+4) .eq. ' zero') then
l=lz+1
line(1:l)='Z'//line(1:lz)
else
line(1:l)=line(lz+1:l)
l=l-lz
endif
endif
if (index(line(i+j:),':') .ne. 0) then
iostat=1
else
read(line(i+j:), *, iostat=iostat) val
endif
endif
if (iostat .eq. 0) then
call putval(line(1:l), val)
else
call putval(line(1:l)//'='//line(i+j:), 0.0)
endif
if (line(1:l) .eq. 'wavelength') then
call dat_group(2, putval)
endif
goto 1
endif
l=index(line, 'Steps:')
if (l .gt. 0) then
do j=1,mcol
values(j)=0
enddo
read(line(l+6:), *,iostat=iostat) values
do j=1,mcol
if (values(j) .ne. 0) then
xcol0=j+1
goto 19
endif
enddo
endif
l=index(line(21:), ',')
if (l .eq. 0) l=10
xaxis=line(21:20+l-1)
19 continue
read(lun, '(a)', err=99,end=99) line
l=index(line, 'Mode: ')
if (l .ne. 0) then
preset=line(l+6:)
l=index(preset, ',')
if (l .eq. 0) goto 99
preset(l:)=' '
call putval('Preset='//preset, 0.0)
l=index(line, 'Preset')
if (l .eq. 0) goto 99
read(line(l+6:), *, err=99,end=99) ymon
pol=.false.
else ! polarized scan
22 read(lun, '(a)', err=99,end=99) line
if (line(1:17) .eq. 'zero for plotting') goto 22
preset='mn'
ymon=0
pol=.true.
endif
ipol=1
yaxis2=' '
call dat_start_options
i=0
call dat_str_option('x', xaxis)
yaxis='Counts'
call dat_str_option('y', yaxis)
call dat_str_option('y2', yaxis2)
if (yaxis .eq. '1') then
if (pol) then
yaxis='up'
else
yaxis='Monitor1'
endif
elseif (yaxis .eq. '2') then
if (pol) then
yaxis='dn'
else
yaxis='Monitor2'
endif
elseif (yaxis .eq. '3') then
yaxis='Monitor3'
else if (yaxis .eq. 'Counts') then
if (pol) then
yaxis='up'
yaxis2='dn'
ipol=2
endif
endif
read(lun,'(a)',err=99,end=99) line
i=1
line(len(line):len(line))=' '
ncol=0
ycol=0
ycol2=0
pcol=0
col2=' '
31 do while (line(i:i) .eq. ' ')
i=i+1
if (i .gt. len(line)) goto 39
enddo
l=i
do while (line(i:i) .ne. ' ')
i=i+1
enddo
ncol=ncol+1
if (ncol .eq. 2) col2=line(l:i)
if (line(l:i) .eq. yaxis .and. ycol .eq. 0) then
ycol=ncol
elseif (line(l:i) .eq. yaxis2 .and. ycol2 .eq. 0) then
ycol2=ncol
elseif (line(l:i) .eq. preset .or.
1 line(l:i) .eq. 'Monitor1' .and. preset .eq. 'Monitor') then
pcol=ncol
elseif (line(l:i) .eq. xaxis) then
xcol=ncol
elseif (xcol0 .eq. ncol) then
xaxis0 = line(l:i)
endif
goto 31
39 if (ycol .eq. 0) goto 99
if (xcol .eq. 0) then
if (xcol0 .ne. 0) then
xcol = xcol0
xaxis = xaxis0
else
xcol=2
xaxis=col2
endif
endif
call dat_group(1, putval)
call putval('XAxis='//xaxis, 0.0)
call putval('YAxis='//yaxis, 0.0)
if (ycol2 .eq. 0) ipol=1
l=min(mcol,max(xcol,pcol,ycol,ycol2))
40 read(lun,'(a)',end=88,err=88) line
if (line .eq. ' ') goto 40
if (line .eq. 'END-OF-DATA') goto 88
read(line,*,err=99,end=99) (values(j),j=1,l)
if (nread .ge. nmax) goto 29
if (pcol .eq. 0) then
if (ymon .eq. 0) ymon=1.
r=ymon
else
r=values(pcol)
if (r .gt. 0) then
if (ymon .eq. 0) ymon=r
else
if (ymon .eq. 0) ymon=1.
r=ymon
endif
endif
f=ymon/r
if (f .le. 0.0) f=1.0
do i=1,ipol
nread=nread+1
xx(nread)=values(xcol)
if (i .eq. 1) then
y=values(ycol)
else
y=values(ycol2)
endif
if (y .gt. 0) then
s=sqrt(y) ! statistical error of detector
else
s=1
endif
yy(nread)=y*f
ss(nread)=s*f
ww(nread)=r
enddo
goto 40
29 print *,'too many points - truncated'
88 close(lun)
if (ipol .gt. 0) then
call fit_dat_table(1, ipol, (nread+ipol-1)/ipol)
endif
call putval('NP', nread*1.0)
call putval('Monitor', ymon)
return
99 nread=-2
rewind lun
print *,'DAT_SICS: error during read'
call putval('Monitor', 0.0)
return
100 nread=-1
rewind lun
call putval('Monitor', 0.0)
end