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