224 lines
5.3 KiB
Fortran
224 lines
5.3 KiB
Fortran
subroutine fit_export(steparg, typ, file)
|
|
! -----------------------------------------
|
|
|
|
implicit none
|
|
|
|
include 'fit.inc'
|
|
|
|
character typ*(*), file*(*)
|
|
real steparg
|
|
|
|
integer ntry
|
|
parameter (ntry=30)
|
|
|
|
character date*20, line*80, filename*256, typup*32, instr*32
|
|
character sample*80
|
|
real start, step, endr, q0, q1, xv, dq
|
|
real sum, suml, stp, xi, best, prec
|
|
|
|
integer i,l,n,j,dig,k,jj,nbest,i0
|
|
integer lun/2/
|
|
integer iran, iostat
|
|
|
|
if (typ .eq. ' ') then
|
|
write (isyswr,'(/X,A,$)') 'Step size (0: automatic):'
|
|
read(isysrd,'(f40.0)',err=999,end=999) step
|
|
write (isyswr,'(/X,A,$)') 'Output file type:'
|
|
read (isysrd, '(A)',err=999,end=999) typup
|
|
call str_upcase(typup, typup)
|
|
else
|
|
step=steparg
|
|
call str_upcase(typup, typ)
|
|
endif
|
|
if (typup .eq. 'D1A') then
|
|
prec=1.0e-4
|
|
else ! DMC, HRPT, LNSP
|
|
prec=1.0e-3
|
|
endif
|
|
call fit_merge(step)
|
|
start=xval(nxmin)
|
|
endr=xval(nxmax)
|
|
if (nxmax .le. nxmin .or. endr .eq. start) then
|
|
write(isyswr,*) 'not enough points to save'
|
|
goto 99
|
|
endif
|
|
if (step .eq. 0) then ! find best step size
|
|
nbest=0
|
|
best=1e30
|
|
iran=12345 ! make random numbers reproducible
|
|
k=nxmax-nxmin
|
|
5 do n=k,2*k ! loop over possible values n
|
|
stp=(endr-start)/n ! step sizes checked between average step size and one half of average
|
|
if (k .gt. ntry*2) then ! try statistically
|
|
sum=0
|
|
suml=ntry*0.15
|
|
do j=1,ntry
|
|
iran=mod(iran*3125,524287) ! quick and dirty random number generator
|
|
i=mod(iran,k)+nxmin+1
|
|
xi=(xval(i)-start)/stp
|
|
sum=sum+abs(xi-nint(xi))
|
|
if (sum .gt. suml) then
|
|
sum=sum*ntry/j
|
|
goto 6
|
|
endif
|
|
enddo
|
|
endif
|
|
sum=0
|
|
suml=(nxmax-nxmin)*0.1
|
|
do i=nxmin+1,nxmax
|
|
xi=(xval(i)-start)/stp
|
|
sum=sum+abs(xi-nint(xi))
|
|
if (sum .gt. suml) then
|
|
sum=sum*k/(i-nxmin)
|
|
goto 6
|
|
endif
|
|
enddo
|
|
sum=sum/suml
|
|
if (sum .lt. 1.0) then
|
|
nbest=n
|
|
goto 7
|
|
endif
|
|
6 if (sum .lt. best) then
|
|
best=sum
|
|
nbest=n
|
|
endif
|
|
|
|
continue
|
|
enddo
|
|
|
|
7 continue
|
|
|
|
if (nbest .eq. 0) stop 'error in FIT_EXPORT'
|
|
step=(endr-start)/nbest
|
|
endif
|
|
if (step .eq. 0) step=1
|
|
|
|
if (file .eq. ' ') then
|
|
write (isyswr,'(/X,A,$)') 'Output file name:'
|
|
read (isysrd, '(A)',err=999,end=999) filename
|
|
else
|
|
filename=file
|
|
endif
|
|
|
|
step=nint(step/prec)*prec ! correct step for output precision
|
|
print *,'step used: ',step
|
|
|
|
i=nxmin
|
|
n=nint((endr-start)/step)+1
|
|
if (npkt+n .gt. maxdat) then
|
|
write(isyswr,*) 'not enough memory to save all points'
|
|
n=maxdat-npkt
|
|
endif
|
|
i=nxmin+1
|
|
k=1
|
|
jj=npkt+1
|
|
if (endr .lt. 0) then
|
|
k=-1
|
|
jj=npkt+n
|
|
endif
|
|
do j=0,n-1
|
|
xv=start+j*step
|
|
do while (xval(i) .lt. xv .and. i .lt. nxmax)
|
|
i=i+1
|
|
enddo
|
|
c cosmetics for x-axis rounding errors
|
|
q0=(xval(i)-xv)/(xval(i)-xval(i-1))
|
|
dq=abs(xv*5E-7/(xval(i)-xval(i-1)))
|
|
if (q0 .lt. dq) then
|
|
q0=0
|
|
elseif (q0 .gt. 1-dq) then
|
|
q0=1
|
|
endif
|
|
q1=1-q0
|
|
YVAL(jj)=YVAL(i-1)*q0+YVAL(i)*q1
|
|
call cvt_real_str(line(2:8), i0, YVAL(jj), 7,0,0,0)
|
|
sig(jj)=sig(i-1)*q0+sig(i)*q1
|
|
jj=jj+k
|
|
enddo
|
|
if (k .lt. 0) then
|
|
start=-xval(nxmax)
|
|
endif
|
|
endr=start+(n-1)*step
|
|
|
|
if (typup .eq. 'DMC' .or. typup .eq. 'HRPT'
|
|
1 .or. typup .eq. 'LNSP') then
|
|
call sys_parse(filename, l, filename, '.dat', 0)
|
|
|
|
call sys_open(lun, filename, 'w', iostat)
|
|
if (iostat .ne. 0) then
|
|
print *,'Can not open ',filename(1:l)
|
|
return
|
|
endif
|
|
date=' '
|
|
call sym_get_str('Date', l, date)
|
|
call sym_get_str('Instrument', jj, instr)
|
|
if (instr(1:jj) .ne. 'DMC' .and. instr(1:jj) .ne. 'HRPT')
|
|
1 call str_trim(instr,typup,jj)
|
|
call sym_get_str('Title', k, line)
|
|
write(lun,'(3a)') instr(1:jj),', ',line(1:k)
|
|
|
|
write (lun,'(a,f9.5,a,f8.3,a,f7.3,3a)')
|
|
1 'lambda=',wavlen,', T=',temp,', dT=',dtemp
|
|
1 ,', Date=''',date(1:l),''''
|
|
|
|
line(1:1)=' '
|
|
call cvt_real_str(line(2:), l, ymon, 8, 0, 6, 1)
|
|
l=l+1
|
|
|
|
call sym_get_str('sample', j, sample)
|
|
if (sample(1:j) .ne. ' ') then
|
|
call str_trim(line
|
|
1 , line(1:l)//', sample="'//sample(1:j)//'"', l)
|
|
endif
|
|
write(lun,'(3f8.3,a)') start, step, endr, line(1:l)
|
|
dig=1
|
|
elseif (typup .eq. 'D1A') then
|
|
call sys_parse(filename, l, filename, '.d1a', 0)
|
|
|
|
call sys_open(lun, filename, 'w', iostat)
|
|
if (iostat .ne. 0) then
|
|
print *,'Can not open ',filename(1:l)
|
|
return
|
|
endif
|
|
|
|
write(lun,'(a,a68)') 'D1A5 Title: ',itit
|
|
write(lun,*)
|
|
write(lun,'(a)') itit
|
|
write(lun,44) n, temp, 0.0, ymon, 0.0
|
|
44 format(i6,f11.3,f10.3,' 1',2f10.1)
|
|
write(lun,'(3f10.4)') start, step, endr
|
|
dig=2
|
|
else
|
|
write(isyswr,*) 'data type ',typup, ' not yet implemented'
|
|
goto 99
|
|
endif
|
|
|
|
l=0
|
|
line=' '
|
|
do j=npkt+1,npkt+n
|
|
call cvt_real_str(line(l+2:l+8), i, YVAL(j), 7,dig-1,0,0)
|
|
l=l+8
|
|
if (l .eq. 80) then
|
|
write(lun,'(a)') line(1:80)
|
|
l=0
|
|
endif
|
|
enddo
|
|
if (l .gt. 0) write(lun,'(a)') line(1:l)
|
|
l=0
|
|
do j=npkt+1,npkt+n
|
|
call cvt_real_str(line(l+2:l+8), i, sig(j), 7,dig,0,0)
|
|
l=l+8
|
|
if (l .eq. 80) then
|
|
write(lun,'(a)') line(1:80)
|
|
l=0
|
|
endif
|
|
enddo
|
|
if (l .gt. 0) write(lun,'(a)') line(1:l)
|
|
call str_trim(fillis, fillis, l)
|
|
write(lun, '(x,3a)') 'Filelist=''',fillis(1:l),''''
|
|
call sym_list(lun, 0, 2
|
|
1,' file monitor instrument title date lambda '
|
|
1//'temp dtemp sample ')
|
|
99 close(lun)
|
|
999 end
|