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

344 lines
7.6 KiB
Fortran

subroutine dat_spec
! -------------------
external dat_spec_desc
external dat_spec_opts
external dat_spec_read
integer dtype/0/
call dat_init_desc(dtype, dat_spec_desc)
call dat_init_opts(dtype, dat_spec_opts)
call dat_init_read(dtype, dat_spec_read)
end
subroutine dat_spec_desc(text)
! ------------------------------
character*(*) text ! (out) description
! type description
! ----------------------------------
text='SPEC spec data format (esrf)'
end
subroutine dat_spec_opts
! ------------------------
print '(x,a)'
1,'from: first dataset (default: 1)'
1,'to: last dataset (default: from)'
1,'x,y,mon: columns to be read (as number or name)'
1,'space: spaces between header items (sls:1,esrf:2)'
1,' '
1,'err: how to calculate error:'
1,' err=s for sqrt(y), this is the default'
1,' err=c for constant value'
1,' err=p for a factor to be multiplied with y'
1,'val: value for err=c and err=p'
end
subroutine dat_spec_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 maxcol, n_nn, n_nl
parameter (maxcol=64,n_nn=1024,n_nl=64)
real values(maxcol)
integer i, i1, i2, idx, l, m, nread0, iset, spc
integer nc, ncol, mcol, xcol, ycol, ycol0
real r, ymon, errvalue
integer iostat
character line*1024, title*132, errtype*1
character xaxis*16, yaxis*16, monam*16, xup*16, yup*16, mup*16
character unam*16, yaxis0*16
character prefix(n_nl)*4, names(n_nn)*16
integer nidx(n_nl)
integer nn, nl, ni, j
read(lun, '(a)',err=100,end=100) line
if (line(1:3) .ne. '#F') goto 100
read(lun, '(a)',err=100,end=100) line
if (line(1:3) .ne. '#E') goto 100
read(lun, '(a)',err=100,end=100) line
if (line(1:3) .ne. '#D') goto 100
read(lun, '(a)',err=100,end=100) title
if (title(1:3) .ne. '#C') goto 100
call putval('Date='//line(4:), 0.0)
call putval('Title='//title(4:), 0.0)
!----- options ------
call dat_start_options
i1=0
call dat_int_option('from', i1)
i2=0
call dat_int_option('to', i2)
spc=0
call dat_int_option('space', spc)
if (spc .le. 0) then
spc=-1
else
spc=spc-1
endif
if (i2 .eq. 0) then
if (i1 .eq. 0) then
i1=1
i2=999999
else
i2=i1
endif
endif
call dat_get_index(idx)
if (idx .ne. 0) then
i1=i1+idx-1
i2=i1
endif
xaxis=' '
call dat_str_option('x', xaxis)
yaxis=' '
call dat_str_option('y', yaxis)
monam='Monitor'
call dat_str_option('mon', monam)
call str_upcase(xup, xaxis)
call str_upcase(yup, yaxis)
call str_upcase(mup, monam)
errtype='s'
call dat_str_option('err', errtype)
call str_upcase(errtype, errtype)
errvalue=1.0
call dat_real_option('val', errvalue)
if (errvalue .le. 0.0) then
print *,'value for error must be > 0'
errvalue=1.0
endif
ymon=0
!----- end options ------
nn=0
nl=0
!--- read parameter names ---
5 read(lun, '(a)', err=99, end=50) line
if (nl .ge. n_nl) goto 9
if (line(1:3) .eq. '#UE') then ! this is not present at ESRF
if (spc .lt. 0) spc = 0 ! assume SLS format (one space as separator)
nl=nl+1
prefix(nl)='UH'//line(4:5)
else if (line(1:2) .eq. '#O') then
nl=nl+1
prefix(nl)='P'//line(3:4)
else if (line(1:1) .eq. '#') then
goto 5
else
goto 9
endif
if (spc .lt. 0) spc=1 ! of no #UE lines, ESRF format (2 spaces as sep.)
nidx(nl)=nn
line(len(line)-spc:)=' ' ! stopper at end
i=5
6 continue
do while (line(i:i) .eq. ' ')
i=i+1
if (i .gt. len(line)) goto 5
enddo
l=i
do while (line(i:i+spc) .ne. ' ')
i=i+1
enddo
if (nn .lt. n_nn) then
nn=nn+1
names(nn)=line(l:i)
endif
goto 6
9 continue
nidx(nl+1)=nn
m=0
nread=0
10 read(lun, '(a)', err=99, end=50) line
if (line(1:3) .ne. '#S') goto 10
read(line(4:), *, err=99, end=99) iset
if (iset .gt. i2) goto 50
if (iset .lt. i1) goto 10
call dat_group(2, putval)
12 read(lun, '(a)', err=99, end=99) line
if (line(1:2) .eq. '#P' .or. line(1:3) .eq. '#UH') then
if (line(1:2) .eq. '#P') then
l=4
else
l=5
endif
do i=1,nl
if (prefix(i) .eq. line(2:l)) then
ni=nidx(i+1)-nidx(i)
do j=1,ni
values(j)=0
enddo
read(line(l:), *, iostat=iostat) (values(j),j=1,ni)
do j=1,ni
call putval(names(nidx(i)+j), values(j))
enddo
endif
enddo
endif
if (line(1:3) .ne. '#N') goto 12
call dat_group(1, putval)
read(line(4:), *, err=99, end=99) ncol
read(lun, '(a)', err=99, end=99) line
if (line(1:3) .ne. '#L') goto 99
xcol=0
mcol=0
ycol=0
ycol0=0
yaxis0=' '
i=3
line(len(line)-1:)=' '
do nc=1,ncol
if (line(i:) .eq. ' ') goto 39
do while (line(i:i) .eq. ' ')
i=i+1
enddo
l=i
31 do while (line(i:i+spc) .ne. ' ')
i=i+1
enddo
call str_upcase(unam, line(l:i-1))
if (unam .eq. xup) then
xcol=nc
xaxis=line(l:i-1)
endif
if (unam .eq. yup) then
ycol=nc
yaxis=line(l:i-1)
endif
if (ycol .eq. 0 .and. ycol0 .eq. 0) then
if (unam .eq. 'DETECTOR' .or. unam .eq. 'APD') then
ycol0=nc
yaxis0=line(l:i-1)
endif
endif
if (unam .eq. mup) mcol=nc
enddo
39 continue
if (xcol .eq. 0) then
read(xup, *, iostat=iostat) xcol
endif
if (ycol .eq. 0) then
read(yup, *, iostat=iostat) ycol
endif
if (mcol .eq. 0) then
read(mup, *, iostat=iostat) mcol
endif
if (xcol .eq. 0) then
xcol=1
if (xaxis .ne. ' ') then
call str_trim(xaxis, xaxis, l)
print *,'DAT_SPEC: ',xaxis(1:l),' not found, take 1st column'
endif
endif
if (ycol .eq. 0) then
if (ycol0 .eq. 0) then
print *,'DAT_SPEC: column not found: ',yaxis
goto 99
endif
if (yup .ne. ' ') then
call str_trim(yaxis, yaxis, l)
print *,'DAT_SPEC: ',yaxis(1:l),' not found, take ',yaxis0
endif
ycol=ycol0
yaxis=yaxis0
endif
call putval('XAxis='//xaxis,0.0)
call putval('YAxis='//yaxis,0.0)
l=min(maxcol,max(mcol,xcol,ycol))
nread0=nread
40 read(lun, *, err=49,end=49) (values(i),i=1,l)
if (nread .ge. nmax) goto 49
nread=nread+1
if (mcol .eq. 0) then
r=1
else
r=values(mcol)
if (ymon .eq. 0) ymon=r
if (r .eq. 0) r=1
endif
ww(nread)=r
yy(nread)=values(ycol)
xx(nread)=values(xcol)
goto 40
49 m=m+1
call fit_dat_table(m, 1, nread-nread0)
goto 10
50 continue
if (errtype .eq. 'S') then
do i=1,nread
ss(i)=sqrt(max(1.0,yy(i)*errvalue))
if (ymon .gt. 0) then
yy(i)=yy(i)*ymon/ww(i)
ss(i)=ss(i)*ymon/ww(i)
endif
enddo
else if (errtype .eq. 'P') then
do i=1,nread
ss(i)=yy(i)*errvalue
if (ss(i) .eq. 0.0) ss(i)=1.0
if (ymon .gt. 0) then
yy(i)=yy(i)*ymon/ww(i)
ss(i)=ss(i)*ymon/ww(i)
endif
enddo
else
if (errtype .ne. 'C') then
print *,'illegal option: err=',errtype
endif
do i=1,nread
ss(i)=errvalue
if (ymon .gt. 0) then
yy(i)=yy(i)*ymon/ww(i)
endif
enddo
endif
call putval('Monitor', ymon)
close(lun)
return
99 print *,'DAT_SPEC: error during read'
98 nread=-2
rewind lun
return
100 nread=-1
rewind lun
end