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

291 lines
6.7 KiB
Fortran

subroutine dat_lnsp
c -------------------
external dat_lnsp_desc
external dat_lnsp_opts
external dat_lnsp_read
integer dtype/0/
call dat_init_desc(dtype, dat_lnsp_desc)
call dat_init_opts(dtype, dat_lnsp_opts)
call dat_init_read(dtype, dat_lnsp_read)
end
subroutine dat_lnsp_desc(text)
! ------------------------------
character*(*) text ! (out) description
! type description
! ----------------------------------
text='LNSP LNS powder Ascii format (DMC, HRPT)'
end
subroutine dat_lnsp_opts
! ------------------------
print '(x,a)'
1,'x: xaxis (2theta,d,Q), default= 2theta'
1,'lambda: used for transformation, if not given if file'
end
subroutine dat_lnsp_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
character line*132, title*80, short*8
character date*20, preset*4, instr*8
character xaxis*64
integer npk
integer i, l, j
real ymon
real meastime, lambda, temp, dtemp
real uur
real xmin, xstep, xend, xval, fact
real yint(0:99)
logical overflow
external dat_lnsp_val
nread=0
preset=' '
overflow=.false.
xstep=0
ymon=0
read(lun,'(a)',err=900,end=900) line
rewind lun
i=index(line, ',')
if (i .gt. 1 .and. i .le. 9) then
instr=line(1:i-1)
else
instr=' '
endif
if (instr .eq. ' ') then
if (line(32:40) .eq. 'Phase No:') then ! fullprof sub output
read(line(1:31), *, err=901,end=901) xmin, xstep, xend
title=line(32:)
read(lun,*,err=901,end=901)
instr='fp'
else
if (forced .le. 0 .and. line(64:67) .ne. 'YTIM') then
inquire(lun, name=line)
call sys_parse(short, l, line, ' ', 2)
call str_upcase(short(1:4), short(1:4))
if (short(1:4) .ne. '.DMC') goto 900
endif
instr='DMC'
endif
elseif (instr .ne. 'DMC' .and. instr .ne. 'HRPT'
1 .and. instr .ne. 'LNSP') then
goto 900
endif
call dat_start_options
xaxis=' '
call dat_str_option('x', xaxis)
lambda=0.0
call dat_real_option('lambda', lambda)
call dat_group(1, putval)
if (instr .eq. 'fp') goto 108
read(lun,'(a)',err=99,end=99) line
if (line(64:67) .eq. 'YTIM') then ! old header
date=' '
meastime=0.0
read(line, '(21x,a20,27x,f8.0)', err=101,end=101) date, meastime
101 read(lun, '(a)',err=901,end=99) title
if (title(67:69) .eq. 'WL:') then
if (lambda .eq. 0.0) then
read(title(70:80), *, err=102,end=99) lambda
endif
title(67:)=' '
102 continue
endif
temp=0
dtemp=0
uur=0
read(lun, '(3f8.0,8x,f8.0,3x,f5.0,8x,f8.0,5x,f8.0)'
1 , err=103,end=99)
1 xmin, xstep, xend, temp, dtemp, uur, ymon
103 if (temp .ne. 0) then
call putval('Temp', temp)
if (dtemp .ne. 0) call putval('dTemp', dtemp)
endif
if (date .ne. ' ') call putval('Date='//date, 0.0)
call dat_group(2, putval)
if (meastime .ne. 0) call putval('MeasTime', meastime)
if (uur .ne. 0) call putval('MuR', uur)
if (ymon .gt. 0) ymon=ymon*1000.
else ! new header
title=line(6:)
read(lun,'(a)',err=99,end=99) line
call dat_delimiters(',','=','''')
call dat_intprt(line, dat_lnsp_val, putval)
read(lun,'(a)',err=99,end=99) line
read(line, *, err=107,end=107) xmin, xstep, xend, ymon
107 j=index(line, ',')
call dat_delimiters(',','=',' ')
if (j .ne. 0 .and. line(j+1:) .ne. ' ') then
call dat_intprt(line(j+1:), dat_lnsp_val, putval)
endif
endif
108 if (xstep .eq. 0) goto 901
if (ymon .gt. 0) then
preset='MN'
elseif (ymon .lt. 0) then
ymon=-ymon
preset='TI'
else
ymon=0. ! was 1.
endif
call dat_group(3, putval)
if (preset .ne. ' ') call putval('Preset='//preset, 0.0)
if (lambda .eq. 0.0) then
call sym_get_real('lambda', lambda)
else if (lambda .ne. 0.0) then
call putval('lambda', lambda)
endif
call dat_powder_init(lambda, '2theta', xaxis)
if (xstep .eq. 0) goto 901
npk=nint((xend-xmin)/xstep)+1
! read intensity block
do i=0,npk-1,10
read(lun, '(10F8.0)', err=902,end=902) (yint(j),j=0,9)
do j=0,min(9,npk-1-i)
if (nread .ge. nmax) then
if (.not. overflow) then
print *,'DAT_LNSP: Too many datapoints, truncated'
overflow=.true.
endif
else
xval=xmin+nread*xstep
nread=nread+1
call dat_powder_trf(xval, xx(nread), fact)
yy(nread)=yint(j)*fact
ww(nread)=fact
if (yint(j) .gt. 0) then
ss(nread)=sqrt(yint(j))
else
ss(nread)=0.0
endif
endif
enddo
enddo
! read sigma block
yint(0)=0
do i=1,npk,10
read(lun, '(10F8.0)', err=290,end=290) (yint(j),j=0,9)
do j=0,9
if (i+j .ge. nread) goto 202
ss(i+j)=yint(j)
enddo
enddo
202 goto 500
290 if (i .le. 1 .and. yint(0) .eq. 0) then
if (instr .ne. 'fp') then
print *,'--- Old LNSP format, no sigma block'
endif
else
print *,'--- Error in sigma block'
endif
500 do i=1,nread
if (ss(i) .eq. 0.0) then
ss(i)=1.0
else
ss(i)=ss(i)*ww(i)
endif
ww(i)=ymon
enddo
call dat_group(3, putval)
if (xaxis .eq. ' ') then
call putval('XAxis=2-Theta [deg]', 0.0)
else
call putval('XAxis='//xaxis, 0.0)
endif
call putval('YAxis=Intensity', 0.0)
call dat_group(1, putval)
call putval('Monitor', ymon)
call str_trim(title, title, l)
call putval('Title='//title(1:l), 0.0)
call putval('Instrument='//instr, 0.0)
call sym_read(lun, dat_lnsp_val)
990 close(lun)
return
! error messages
900 nread=-1
rewind lun
return
901 print *,'DAT_LNSP: Error in header'
goto 990
902 print *,'DAT_LNSP: Error in intensity block'
goto 990
99 print *,'DAT_LNSP: error during read'
rewind lun
nread=-2
return
end
subroutine dat_lnsp_val(str, val, putval)
character*(*) str
real val
external putval
integer i
if (val .eq. 0) then
i=index(str, '=')
else
i=0
endif
if (i .eq. 0) then ! numeric
if (str .eq. 'T') then
call putval('Temp', val)
return
elseif (str .eq. 'dT') then
call putval('dTemp', val)
return
elseif (str .eq. ' ') then
return
endif
else
if (str(1:i) .eq. 'Filelist=') return
endif
call putval(str, val)
end