445 lines
11 KiB
Fortran
445 lines
11 KiB
Fortran
program polcal
|
|
|
|
c this program work, because the dat_tasmad_read routine is
|
|
c coding the polarisation information into the x-value ip.b
|
|
c where i is 1..9 for xx, yx, yz, zy, yy, zy, xz, yz, zz
|
|
c p is 1..4 for ++, +-, -+, --
|
|
c b is 0 for signal and 1 for background
|
|
|
|
implicit none
|
|
character files*1024
|
|
character bfiles*1024
|
|
character line*2048
|
|
character word*1024
|
|
|
|
integer l, pos
|
|
integer i
|
|
real xx(36), yy(36), ss(36)
|
|
real cnts(9,4,2), sigma(9,4,2)
|
|
integer n
|
|
character mukind*20
|
|
character single*1024
|
|
logical merge,bgr
|
|
integer listflag
|
|
real numor, mon, lastnumor
|
|
integer status
|
|
integer iostat
|
|
integer summary,sumx,sumy
|
|
common summary,sumx,sumy
|
|
|
|
call sys_get_cmdpar(line, l)
|
|
if (line .eq. ' ') then
|
|
print *
|
|
print *,'polcal calculates polarisation matrices from files ',
|
|
1 'measured with the'
|
|
print *,'polmat command on TASP/MuPAD'
|
|
print *
|
|
print *,'Usage:'
|
|
print *,
|
|
1 'polcal <numors> [-m] [-b <background numors>]'
|
|
print *
|
|
print *,'Options:'
|
|
print *
|
|
print *,' -m'
|
|
print *,' merge files'
|
|
print *,' without this option, for every file a polarisation'
|
|
print *,' matrix is calculated'
|
|
print *
|
|
print *,' -b'
|
|
print *,' treat numors given before -b as signal,'
|
|
print *,' numors after -b as background.'
|
|
print *,' the files are merged even without the -m option.'
|
|
print *,' without this option, polcal sorts out automatically'
|
|
print *,' signal and background files'
|
|
print *
|
|
print *,' -f'
|
|
print *,' make summary files fort.11 etc. with the'
|
|
print *,' matrix elements'
|
|
print *
|
|
print *,' -f1213'
|
|
print *,' make a summary of Pxy vs Pxz (output file fort.1)'
|
|
print *
|
|
goto 9
|
|
endif
|
|
|
|
c Argument processing
|
|
pos=1
|
|
files=' '
|
|
merge=.false.
|
|
summary=0
|
|
bgr=.false.
|
|
call str_get_word(line, pos, word)
|
|
1 do while (word .ne. ' ')
|
|
if (word .eq. '-b') then
|
|
call str_get_word(line, pos, bfiles)
|
|
bgr=.true.
|
|
merge=.true.
|
|
elseif (word .eq. '-m') then
|
|
merge=.true.
|
|
elseif (word(1:2) .eq. '-f') then
|
|
if (word(3:) .eq. ' ') then
|
|
summary=1
|
|
else
|
|
summary=2
|
|
sumx=0
|
|
sumy=0
|
|
read(word(3:4), *, iostat=iostat) sumx
|
|
read(word(5:6), *, iostat=iostat) sumy
|
|
endif
|
|
else
|
|
if (files .eq. ' ') then
|
|
files = word
|
|
else
|
|
call str_trim(files, files, l)
|
|
if (l .lt. len(files) - 1) then
|
|
l=l+1
|
|
files(l:l)=','
|
|
files(l+1:) = word
|
|
endif
|
|
endif
|
|
endif
|
|
call str_get_word(line, pos, word)
|
|
enddo
|
|
|
|
call sys_setenv('dat_defspec','TASP')
|
|
call fit_init_silent
|
|
call clr_me(cnts)
|
|
call clr_me(sigma)
|
|
|
|
if (merge) then
|
|
if (bgr) then
|
|
call fit_dat_silent
|
|
call fit_dat(files)
|
|
call fit_merge(0.2) ! merge signal and bgr
|
|
call fit_get_real('Monitor', mon)
|
|
call fit_get_array('X', xx, 36, n)
|
|
call fit_get_array('Y', yy, 36, n)
|
|
call fit_get_array('S', ss, 36, n)
|
|
call fill_me(cnts, xx, yy, n, 1) ! 1 = force to signal
|
|
call fill_me(sigma, xx, ss, n, 1)
|
|
call calc_pol(cnts, sigma, status)
|
|
if (status .eq. 0) then
|
|
print *,'NO DATA'
|
|
goto 9
|
|
else
|
|
call info_line(files, 'S', 1)
|
|
endif
|
|
call fit_dat_silent
|
|
call fit_dat(bfiles)
|
|
call fit_merge(0.2) ! merge signal and bgr
|
|
call fit_mon(mon)
|
|
call fit_get_array('X', xx, 36, n)
|
|
call fit_get_array('Y', yy, 36, n)
|
|
call fit_get_array('S', ss, 36, n)
|
|
call fill_me(cnts, xx, yy, n, 2) ! 2 = force to bgr
|
|
call fill_me(sigma, xx, ss, n, 2)
|
|
call calc_pol(cnts, sigma, status)
|
|
if (status .eq. 0) then
|
|
print *,'NO DATA'
|
|
else
|
|
call info_line(bfiles, 'B', 2)
|
|
endif
|
|
else
|
|
call fit_dat_silent
|
|
call fit_dat(files)
|
|
call fit_merge(0.01) ! do not merge signal and bgr
|
|
call fit_get_array('X', xx, 36, n)
|
|
call fit_get_array('Y', yy, 36, n)
|
|
call fit_get_array('S', ss, 36, n)
|
|
call fill_me(cnts, xx, yy, n, 0)
|
|
call fill_me(sigma, xx, ss, n, 0)
|
|
call calc_pol(cnts, sigma, status)
|
|
if (status .eq. 0) then
|
|
print *,'NO DATA'
|
|
else
|
|
call info_line(files, 'S', 1)
|
|
if (status .eq. 2) call info_line(files, 'B', 1)
|
|
endif
|
|
endif
|
|
else
|
|
listflag=0
|
|
lastnumor = -1
|
|
do i=1,99999
|
|
call fit_dat_silent
|
|
call fit_dat_next(files, listflag)
|
|
if (listflag .eq. 0) goto 9
|
|
mukind=' '
|
|
call fit_get_str('mukind', l, mukind)
|
|
if (mukind .ne. 'background') then
|
|
call clr_me(cnts)
|
|
call clr_me(sigma)
|
|
if (mukind .eq. ' ') goto 19
|
|
endif
|
|
call fit_get_array('X', xx, 36, n)
|
|
call fit_get_array('Y', yy, 36, n)
|
|
call fit_get_array('S', ss, 36, n)
|
|
call fill_me(cnts, xx, yy, n, 0)
|
|
call fill_me(sigma, xx, ss, n, 0)
|
|
call calc_pol(cnts, sigma, status)
|
|
if (status .eq. 0) goto 19
|
|
call fit_get_real('Numor', numor)
|
|
if (mukind .eq. 'background') then
|
|
if (lastnumor .eq. -1) then
|
|
print *,'WARNING: background only'
|
|
call cvtnumor(single, numor)
|
|
call info_line(single, ' ', 1)
|
|
else
|
|
call cvtnumor(single, lastnumor)
|
|
call cvtnumor(single, numor)
|
|
call info_line(single, 'B', 2)
|
|
endif
|
|
else
|
|
call cvtnumor(single, numor)
|
|
call info_line(single, ' ', 1)
|
|
lastnumor=numor
|
|
endif
|
|
19 continue
|
|
enddo
|
|
endif
|
|
9 end
|
|
|
|
|
|
subroutine cvtnumor(result, numor)
|
|
character result*(*)
|
|
real numor
|
|
integer l
|
|
|
|
write(result,'(f10.0)') numor
|
|
call str_first_nonblank(result, l)
|
|
result = result(l:9)
|
|
end
|
|
|
|
|
|
subroutine info_line(file, type, prthold)
|
|
|
|
character file*(*), type*(*)
|
|
integer prthold
|
|
|
|
character line1*132/' '/, line2*132/' '/
|
|
real vmin(5), vmax(5), mean(5), diff(5)
|
|
integer i,l
|
|
character line*132, label*4
|
|
logical secondline
|
|
|
|
if (prthold .eq. 2) then
|
|
if (line1(1:3) .eq. ' ') line1(1:3)='sig'
|
|
call str_trim(line1, line1, l)
|
|
print '(x,a)',line1(1:l)
|
|
call str_trim(line2, line2, l)
|
|
if (l .gt. 1) print '(x,a)',line2(1:l)
|
|
endif
|
|
|
|
do i=1,5
|
|
vmin(i)=0
|
|
vmax(i)=0
|
|
enddo
|
|
if (type .eq. 'B') then
|
|
call meta_real_range('QH', vmin(1), vmax(1))
|
|
call meta_real_range('QH_B', vmin(1), vmax(1))
|
|
call meta_real_range('QK', vmin(2), vmax(2))
|
|
call meta_real_range('QK_B', vmin(2), vmax(2))
|
|
call meta_real_range('QL', vmin(3), vmax(3))
|
|
call meta_real_range('QL_B', vmin(3), vmax(3))
|
|
call meta_real_range('EN', vmin(4), vmax(4))
|
|
call meta_real_range('EN_B', vmin(4), vmax(4))
|
|
call meta_real_range('TEMP', vmin(5), vmax(5))
|
|
call meta_real_range('TEMP_B', vmin(5), vmax(5))
|
|
label='bgr'
|
|
else
|
|
call meta_real_range('QH', vmin(1), vmax(1))
|
|
call meta_real_range('QH_S', vmin(1), vmax(1))
|
|
call meta_real_range('QK', vmin(2), vmax(2))
|
|
call meta_real_range('QK_S', vmin(2), vmax(2))
|
|
call meta_real_range('QL', vmin(3), vmax(3))
|
|
call meta_real_range('QL_S', vmin(3), vmax(3))
|
|
call meta_real_range('EN', vmin(4), vmax(4))
|
|
call meta_real_range('EN_S', vmin(4), vmax(4))
|
|
call meta_real_range('TEMP', vmin(5), vmax(5))
|
|
call meta_real_range('TEMP_S', vmin(5), vmax(5))
|
|
if (type .eq. 'S') then
|
|
label='sig'
|
|
else
|
|
label=' '
|
|
endif
|
|
endif
|
|
|
|
do i=1,5
|
|
mean(i) = 0.5 * (vmin(i) + vmax(i))
|
|
diff(i) = 0.5 * (vmax(i) - vmin(i))
|
|
enddo
|
|
! temperature tolerance 10 %
|
|
secondline = (diff(5) .gt. mean(5) * 0.1)
|
|
! q tolerance 0.003
|
|
do i=1,3
|
|
if (diff(i) .gt. 0.003) secondline = .true.
|
|
enddo
|
|
! en tolerance 0.01
|
|
if (diff(4) .gt. 0.01) secondline = .true.
|
|
|
|
call str_trim(file, file, l)
|
|
write(line1, '(2a,3f8.3,'' en:'',f8.2,'' T:'',f8.3,2a)')
|
|
1 label, 'q:', mean, ' file(s): ',file(1:l)
|
|
|
|
line2=' '
|
|
if (secondline) then
|
|
line=' '
|
|
do i=1,5
|
|
if (i .eq. 4) then
|
|
write(line(i*8-7:i*8), '(f8.2)') diff(i)
|
|
else
|
|
write(line(i*8-7:i*8), '(f8.3)') diff(i)
|
|
endif
|
|
enddo
|
|
write(line2,'(7a)') ' +/-',line(1:24),' +/-',line(25:32)
|
|
1 ,' +/-',line(33:40),' <-- mixed values!!!'
|
|
endif
|
|
|
|
if (prthold .gt. 0) then
|
|
call str_trim(line1, line1, l)
|
|
print '(x,a)',line1(1:l)
|
|
call str_trim(line2, line2, l)
|
|
if (l .gt. 1) print '(x,a)',line2(1:l)
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine calc_pol(cnts, sigma, status)
|
|
implicit none
|
|
real cnts(9,4,2), sigma(9,4,2)
|
|
integer status
|
|
|
|
integer i,i1,i2,k,ip,im
|
|
real sum2, d2, sum, pij, dpij2, c1, c2, sq1, sq2
|
|
character line*80
|
|
|
|
real pm(3,4), sm(3,4)
|
|
logical done,bgr,matoutput
|
|
integer l
|
|
integer summary,sumx,sumy
|
|
common summary,sumx,sumy
|
|
real numor
|
|
|
|
status=0
|
|
matoutput = .false.
|
|
c k=1: normal polarity, k=2: NEG. polarity
|
|
do k=1,2
|
|
done = .false.
|
|
bgr = .false.
|
|
ip=k*2-1
|
|
im=k*2
|
|
do i1=1,3
|
|
sum2 = 0
|
|
d2 = 0
|
|
line=' '
|
|
do i2=1,3
|
|
i = (i1-1)*3 + i2
|
|
if (cnts(i,ip,2) .gt. 0 .or. cnts(i,im,2) .gt. 0) status=2
|
|
c1 = cnts(i, ip, 1) - cnts(i, ip, 2)
|
|
c2 = cnts(i, im, 1) - cnts(i, im, 2)
|
|
sq1 = sigma(i, ip, 1)**2 + sigma(i, ip, 2)**2
|
|
sq2 = sigma(i, im, 1)**2 + sigma(i, im, 2)**2
|
|
sum = c1 + c2
|
|
if (sum .eq. 0 .or.
|
|
1 cnts(i, ip, 1) .eq. 0 .and. cnts(i, im, 1) .eq. 0) then
|
|
pij=0
|
|
dpij2 = 1
|
|
else
|
|
done = .true.
|
|
pij = (c1 - c2) / sum
|
|
dpij2 = ((1-pij)/sum) ** 2 * sq1
|
|
dpij2 = dpij2 + ((1+pij)/sum) ** 2 * sq2
|
|
endif
|
|
pm(i1, i2) = pij
|
|
sm(i1, i2) = sqrt(dpij2)
|
|
sum2 = sum2 + pij*pij
|
|
d2 = d2 + dpij2 * pij * pij
|
|
enddo
|
|
pm(i1,4) = sqrt(sum2)
|
|
if (sum2 .eq. 0) then
|
|
sm(i1,4)=1
|
|
else
|
|
sm(i1,4) = sqrt(d2/sum2)
|
|
endif
|
|
enddo
|
|
if (done) then
|
|
matoutput = .true.
|
|
if (status .eq. 0) status=1
|
|
print *,
|
|
1'------------------------------------------------------'
|
|
call fit_get_str('Title', l, line)
|
|
if (k .eq. 2) then
|
|
print *,'NEG. polarity ',line(1:l)
|
|
else
|
|
print *,'normal polarity ',line(1:l)
|
|
endif
|
|
print *,' pix sigma piy sigma piz sigma'
|
|
1, ' |Pi| sigma'
|
|
do i1=1,3
|
|
line=' '
|
|
if (summary .eq. 1) call fit_get_real('numor', numor)
|
|
do i2=1,4
|
|
if (pm(i1, i2) .ne. 0 .or. sm(i1, i2) .ne. 1.0) then
|
|
write(line(i2*20-19:i2*20), '(2f7.3)') pm(i1, i2), sm(i1, i2)
|
|
endif
|
|
if (summary .eq. 1) then
|
|
write(i1*10+i2,*) numor, pm(i1,i2), sm(i1,i2)
|
|
endif
|
|
enddo
|
|
print *,line(1:78)
|
|
if (summary .eq. 2) then
|
|
write(1,*) pm(sumx/10, mod(sumx,10))
|
|
1 , pm(sumy/10, mod(sumy,10)),sm(sumy/10,mod(sumy,10))
|
|
endif
|
|
enddo
|
|
endif
|
|
enddo
|
|
if (.not. matoutput) then
|
|
print *,'WARNING: zero count matrix can not be calculated'
|
|
endif
|
|
end
|
|
|
|
subroutine clr_me(values)
|
|
real values(9,4,2)
|
|
integer i,j,k
|
|
do k=1,2
|
|
do j=1,4
|
|
do i=1,9
|
|
values(i,j,k)=0
|
|
enddo
|
|
enddo
|
|
enddo
|
|
end
|
|
|
|
subroutine fill_me(values, xx, yy, n, mode)
|
|
|
|
implicit none
|
|
integer n
|
|
integer mode
|
|
real values(9,4,2), xx(n), yy(n)
|
|
|
|
integer i,j,k,m,dbl
|
|
|
|
dbl=0
|
|
do m=1,n
|
|
i = nint(xx(m))
|
|
if (mode .ne. 0) then
|
|
k = mode
|
|
else if (abs(xx(m)-i) .lt. 0.05) then
|
|
k = 1
|
|
else
|
|
k = 2
|
|
endif
|
|
j=mod(i,10)
|
|
i=i/10
|
|
if (i .gt. 0 .and. i .le. 9 .and. j .gt. 0 .and. j .le. 4) then
|
|
if (values(i,j,k) .ne. 0) then
|
|
dbl=1
|
|
endif
|
|
values(i,j,k) = yy(m)
|
|
endif
|
|
enddo
|
|
if (dbl .ne. 0) then
|
|
print *,'WARNING: duplicate points skipped'
|
|
endif
|
|
end
|