	program epsilon
	COMMon/namen/fnaam
	character*80 fnaam
	open(34,file='calibratie.txt')
	
	open(33,file='eps_study.txt')
	do i=24,24
	  if (i.eq.0) fnaam='../0808/cali_eps_0808.txt'
	  if (i.eq.1) fnaam='../bol1/cali_eps_bol1.txt'
	  if (i.eq.2) fnaam='../bol2/cali_eps.txt'
	  if (i.eq.3) fnaam='../bol3/cali_eps.txt'
	  if (i.eq.4) fnaam='../bol4/cali_eps.txt'
	  if (i.eq.5) fnaam='../bol5/cali_eps.txt'
	  if (i.eq.6) fnaam='../bol6/cali_eps.txt'
	  if (i.eq.7) fnaam='../bol7/cali_eps.txt'
	  if (i.eq.8) fnaam='../bol8/cali_eps.txt'
	  if (i.eq.9) fnaam='../bol9/cali_eps.txt'
	  if (i.eq.10) fnaam='../bol10/cali_eps.txt'
	  if (i.eq.11) fnaam='../bol11/cali_eps.txt'
	  if (i.eq.12) fnaam='../bol12/cali_eps.txt'
	  if (i.eq.13) fnaam='../bol13/cali_eps.txt'
	  if (i.eq.14) fnaam='../bol14/cali_eps.txt'
	  if (i.eq.15) fnaam='../bol15/cali_eps.txt'
	  if (i.eq.16) fnaam='../bol16/cali_eps.txt'
	  if (i.eq.17) fnaam='../bol17/cali_eps.txt'
	  if (i.eq.18) fnaam='../bol18/cali_eps.txt'
	  if (i.eq.19) fnaam='../bol19/cali_eps.txt'
	  if (i.eq.24) fnaam='../bol24/cali_eps.txt'
	  write(33,*)fnaam
          write(34,*)'CALIBRATIE CONTROLE VOOR ',fnaam
	  call go
	enddo
	close(33)
	close(34)

	end


	subroutine go
	include 'Online.cmm'
	include 'getcor.cmm'
	common /kang/anw,ane,azw,aze,dzdy
        COMMon/messa/imes,messages
        character*100 messages(200)
        real rotz_null(100)
        real dzdy_null(100)
	real dzdy_save(100)
	real dzdy(100)
        real dy_twee(2,8)
	idatum=20000808 ! a bol2 date

	call readeps

	do i=1,100
	   rotz_null(i)=9999999.
	   dzdy_null(i)=9999999.
	enddo
C--
C From Bol0 ....
C--
	rotz_null(1)=-4.8
	rotz_null(2)=17.37
	rotz_null(3)=6.28
	rotz_null(4)=28.96
	rotz_null(5)=-4.8
	rotz_null(6)=17.37
	rotz_null(7)=6.28
	rotz_null(8)=28.96

        dzdy_null(1)=0.0461/1000.
        dzdy_null(2)=0.2345/1000.
        dzdy_null(3)=0.0807/1000.
        dzdy_null(4)=0.2576/1000.
        dzdy_null(5)=0.0154/1000.
        dzdy_null(6)=0.2230/1000.
        dzdy_null(7)=0.0308/1000.
        dzdy_null(8)=0.1884/1000.
C--
C From BOL2 = eerste kamer waar ik vertrouwen in heb
C--
      dy_twee(1,1)=   0.001
      dy_twee(1,2)=  0.002 
      dy_twee(1,3)= -0.003 
      dy_twee(1,4)=  0.003 
      dy_twee(1,5)=  0.001 
      dy_twee(1,6)=  0.004 
      dy_twee(1,7)= -0.001 
      dy_twee(1,8)=  0.003 
      dy_twee(2,1)= 0.007
      dy_twee(2,2)=0.002
      dy_twee(2,3)=0.006
      dy_twee(2,4)=0.004
      dy_twee(2,5)=0.007
      dy_twee(2,6)=0.002
      dy_twee(2,7)=0.003
      dy_twee(2,8)=0.002


C--
C
C--
	
	print *,'*********************************************'
	print *,'DY simple: assume that <dz>==0 --> angles '
	print *,'*********************************************'
	do id=1,8 !  A+ study
	   dzdy(id)=0.
	enddo

	do id=1,8 !  A+ study
	   dzdy(id)=(dzeps(id,2)-dzeps(id,1))/26.011
	   dzdy(id)=dzdy(id)+(dzeps(id,3)-dzeps(id,2))/26.011
	   dzdy(id)=dzdy(id)+2.*(dzeps(id,3)-dzeps(id,1))/(26.011*2.)
           dzdy(id)=dzdy(id)/4.
	enddo

       write(99,*)dzeps(1,2),dzeps(5,2),dzeps(7,2)

	print *,'NW NE SW SE'
	do id=1,8
	   step12=(dyeps(id,2)-dyeps(id,1))-15.017*dzdy(id) ! 
	   step23=(dyeps(id,3)-dyeps(id,2))+15.017*dzdy(id) ! 
	   step13=dyeps(id,3)-dyeps(id,1)
           if (abs(step12).gt.0.010 .and. 
     +             abs(step12-dy_twee(1,id)).gt.0.010) then
	      write(34,555)id,step12
 555       format('Y STEP 1-2 te groot id=',I2,' dy=',F7.3)
           endif
           if (abs(step23).gt.0.010 .and.
     +             abs(step23-dy_twee(2,id)).gt.0.010) then
              write(34,556)id,step23
 556       format('Y STEP 2-3 te groot id=',I2,' dy=',F7.3)
           endif
           if (abs(step13).gt.0.010 .and.
     +             abs(step13-
     +             dy_twee(1,id)-dy_twee(2,id)).gt.0.010) then
              write(34,557)id,step13
 557       format('Y STEP 1-3 te groot id=',I2,' dy=',F7.3)
           endif




	   alf1=(dzeps(id,2)-dzeps(id,1))/26.011
	   alf2=(dzeps(id,3)-dzeps(id,2))/26.011
	   alf3=(dzeps(id,3)-dzeps(id,1))/(26.011*2.)
	   write(*,1000)id,step12,step23,step13,
     +     alf1*1000,alf2*1000,alf3*1000,dzdy(id)*1000
 1000	   format('id',1X,I2,' step 12 23 13 ',3F7.3,' angles: ',4F7.4)


           if (alf1*1000.gt.2 .or. alf2*1000.gt.2 .or.
     +           alf3*1000.gt.2) then
              write(34,558)id,alf1*1000,alf2*1000,alf3*1000
 558       format('dZ hoeken te groot id=',I2,' hoeken=',3F7.3)
           endif

	enddo

C Warnings for inconsistent calibration, angle step from A+, A-
	print *,'NW NE SW SE'
	do id=1,4
	   if (id.eq.1) idp=7
	   if (id.eq.2) idp=8
	   if (id.eq.3) idp=5
	   if (id.eq.4) idp=6
           alf3p=(dzeps(id,3)-dzeps(id,1))/(26.011*2.)
           alf3m=(dzeps(idp,3)-dzeps(idp,1))/(26.011*2.)
           dalf3=alf3p+alf3m
           if (id.eq.2 .or. id.eq.4) dalf3=dalf3-0.4/1000.
	   if (abs(dalf3*1000).gt.0.5) then
              write(34,559)id,idp,dalf3*1000
 559       format('dZ hoek stap A+/A- te groot ',
     +                   I2,I2,' stap=',F7.3)
           endif
         enddo


C store these dzdy, because they get overwritten in next routine
	do i=1,100
		dzdy_save(i)=dzdy(i)
	enddo

	print *,'*********************************************'
	print *,' '
	print *,' '
	print *,'*********************************************'
	print *,'DZ and DY from assumption dz_nw_13 is know' 
	print *,'*********************************************'
	call deltaz_final
	print *,'*********************************************'
	print *,' '
	print *,'*********************************************'
	print *,'Absolute epsilons in A layer 2'
	print *,'Ignore effect of  y=z/angle  and z=-y/angle '
	do id=1,8
	     write(*,3010)id,dzeps(id,2),dyeps(id,2)
	     write(33,3010)id,dzeps(id,2),dyeps(id,2)
 3010      format('dzeps   dyeps ',I2,1X,2F7.3)
        if (abs(dyeps(id,2)).gt.0.500) then 
            write(34,*)'xplate Y verdacht op ',id,dyeps(id,2)
        endif
	enddo 

	do id=1,4
          if (id.eq.1) idp=7
          if (id.eq.2) idp=8
          if (id.eq.3) idp=5
          if (id.eq.4) idp=6
          dy2=dyeps(id,2)+dyeps(idp,2)
	  if (abs(dy2).gt.0.050) then
            write(34,*)'dY stap A+ A- in L2 te groot ',id,idp,dy2
          endif
         enddo
		

        print *,' '
        print *,'*********************************************'
        print *,'Corrected epsilons in A layer 2 '
        do id=1,8
	do il=2,2
	   isi=1
	   ycorp=100.*(dzdy_save(id)-dzdy_null(id))
	   if (id.eq.1 .or. id.eq.2) ycorp=-ycorp
	   if (id.eq.5 .or. id.eq.6) ycorp=-ycorp
           if (id.eq.2 .or. id.eq.3) isi=-1
           if (id.eq.6 .or. id.eq.7) isi=-1
	   if (id.gt.8) print *,'id not implemented ',id
	     ycor=100.*(rotzplot(id)-rotz_null(id))/1000.
             write(*,3011)id,dzeps(id,il),dyeps(id,il)+isi*ycor,
     +                   dyeps(id,il)+isi*ycorp
             write(33,3011)id,dzeps(id,il),dyeps(id,il)+isi*ycor,
     +              dyeps(id,il)+isi*ycorp
 3011      format('dzeps  dyeps_rotz dyeps_ang',I2,3(1X,F7.3))
        if (id.eq.1 .or. id.eq.3 .or.
     +      id.eq.5 .or. id.eq.7) then
	   if (abs(dzeps(id,il)).gt.0.050) then
	    write(34,*)'Ref: Absolute Epsilon Z te groot ',
     +                   id,dzeps	
           endif
        else ! non ref
	  if (abs(dzeps(id,il)).gt.1.000) then
	    write(34,*)'Zeer lange Xplate?  ',id,dzeps	
           endif
        endif
        if (abs(dyeps(id,il)+isi*ycor).gt.0.050 .and. 
     +          abs(dyeps(id,il)+isi*ycorp).gt.0.050) then
	    write(34,*)'Absolute Epsilon Y in L2 to large ',
     +             id,(dyeps(id,il)+isi*ycorp)
        endif
        enddo
        enddo

	end

	subroutine deltaz_fancy(il)
        include 'Online.cmm'
        include 'getcor.cmm'
	common /kang/anw,ane,azw,aze,dzdy
	real dzdy(100)
	if (il.gt.2) stop
        dy=26.011
	s=-1.
        dnw=0.006 ! assumed
C	if (il.eq.2) then
C	  dnw=dzeps(1,il+1)-dzeps(1,il)-dy*anw
C	endif
        anw_0=(dzeps(1,il+1)-dzeps(1,il) - dnw)/dy
	azw=s*(dzeps(5,il+1)-dzeps(5,il) - dnw)/dy
	dzw=dzeps(3,il+1)-dzeps(3,il) - azw*dy
 	anw_1=s*(dzeps(7,il+1)-dzeps(7,il) - dzw)/dy
	anw=(anw_0+anw_1)/2.
	dzw_1=dzeps(7,il+1)-dzeps(7,il) - s*anw*dy
c	dzw=(dzw+dzw_1)/2.
	ane=(dzeps(2,il+1)-dzeps(2,il) - dnw)/dy
	aze=(dzeps(4,il+1)-dzeps(4,il) - dzw)/dy
	dnw_1=(dzeps(6,il+1)-dzeps(6,il) - s*aze*dy)
	aze_1=s*(dzeps(6,il+1)-dzeps(6,il) - dnw)/dy
        aze=(aze+aze_1)/2.
	dzw_2=(dzeps(8,il+1)-dzeps(8,il) - s*ane*dy)
	ane_1=s*(dzeps(8,il+1)-dzeps(8,il) - dzw)/dy
	ane=(ane+ane_1)/2.

	istep=12
	if (il.eq.2) istep=23
	print *,'delta z for dnw=',dnw,' and step ',istep
	write(*,3000)dnw,dzw
	write(*,3001)dnw_1,dzw_1,dzw_2
	write(*,3002)anw,ane,azw,aze
	write(*,3003)anw_1,ane_1,aze_1
 3000	format(' dz NW SW ',2F7.4)
 3001	format(' alternative dz for NW SW SW ',3F7.4)
 3002	format(' angle NW NE SW SE ',4F7.4)
 3003	format(' alternative angle NW NE SE ',4F7.4)
	
       print *,'dy  NW NE SW SE with these angles'
	dzdy(1)=anw
	dzdy(2)=ane
	dzdy(3)=azw
	dzdy(4)=aze
	dzdy(5)=s*azw	
	dzdy(6)=s*aze	
	dzdy(7)=s*anw	
	dzdy(8)=s*ane	
        do id=1,8
           step12=(dyeps(id,2)-dyeps(id,1))-15.017*dzdy(id) !
           step23=(dyeps(id,3)-dyeps(id,2))+15.017*dzdy(id) !
           step13=dyeps(id,3)-dyeps(id,1)
C          print *,'id ',id,' alfa ',dzdy(id)
           write(*,5000)id,step12,step23,dzdy(id)*1000
 5000   format('id',1X,I2,' step 12 23 13 ',2F7.3,' used angle: ',F7.4)
        enddo



	end
	subroutine deltaz_assume(idr,idn)
        include 'Online.cmm'
        include 'getcor.cmm'

* set step 1 to 3 on ref side Z-eps to zero
	alf_ref=(dzeps(idr,3)-dzeps(idr,1))/26./2. ! delta_z:=0

* delta12 ref
	dv=dzeps(idr,2)-dzeps(idr,1)
        delta12_ref=dv-26*alf_ref

* delta23 ref
	dv=dzeps(idr,3)-dzeps(idr,2)
	delta23_ref=dv-26*alf_ref

* alf_non
	dv=dzeps(idn,3)-dzeps(idn,2)
	alf23_non=(dv-delta23_ref)/26.
	dv=dzeps(idn,2)-dzeps(idn,1)
	alf12_non=(dv-delta12_ref)/26.

	print *,'alf_ref -assume ',alf_ref
	print *,'alf_non 23',alf23_non
	print *,'alf_non 12',alf12_non
	print *,'delta12 ',delta12_ref
	print *,'delta23 ',delta23_ref


	end

	subroutine afewpars(idr,idn)
* attempt to calculate epsilons from linear equation: failed, stuff is singlular* as a posteriori expected.
	integer IR(100)
	include 'Online.cmm'
        include 'getcor.cmm'
	integer ilh,irh
	parameter(irh=6)
	parameter(ilh=4)
	double precision v(ilh,irh),dv(irh),av(irh),vin(irh,ilh)
	ifree=4
	do i=1,ilh
	  do j=1,irh
	    dv(j)=0.
	    av(j)=0.
	    v(i,j)=0.d0
	    vin(j,i)=0.d0
	  enddo
	enddo
	
	print *,'hallo 1' 
	vin(1,1)=1.
	vin(1,2)=26. ! should be 26, but then matrix is singular-> rotzRef:=0

        vin(2,1)=1.
	vin(2,2)=0.
	vin(2,3)=26.
	vin(2,4)=0.

	vin(3,2)=26.
	vin(3,4)=1.

        vin(4,3)=26.
	vin(4,4)=1.

	vin(5,1)=1.
	vin(5,2)=52.
	vin(5,3)=0.
	vin(5,4)=1.

	vin(6,1)=1.

	print *,'hallo 2' 
  	do i=1,ilh
  	   do j=1,irh
  		v(i,j)=vin(j,i)
  	   enddo
          enddo

	print *,'hallo 3' 



	dv(1)=dzeps(idr,2)-dzeps(idr,1)
	dv(2)=dzeps(idn,2)-dzeps(idn,1)
	dv(3)=dzeps(idr,3)-dzeps(idr,2)
	dv(4)=dzeps(idn,3)-dzeps(idn,2)
	dv(5)=(dzeps(idr,3)-dzeps(idr,1))
	dv(6)=0.

	print *,'for idr and idn: ',idr,idn
	do i=1,4
	  print *,'dzeps i ',i , dv(i)
        enddo 
	print *,'EQUATION '
        call DSEQN(ilh,v,irh,IR,IFAIL,1,dv)

	print *,' IFAIL ', IFAIL
	print *,' d12 ',dv(1)
	print *,' a_r ',dv(2)
	print *,' a_n ',dv(3)
	print *,' d23 ',dv(4)

	end




	
	subroutine readeps
	COMMon/namen/fnaam
	include 'Online.cmm'
C-- read calibration files
        include 'getcor.cmm'
        character*80 fnaam
C-- some initial stuff
        do i=1,100
           do j=1,3
	      dyeps_save(i,j)=0.
	      dzeps_save(i,j)=0.
	      counts_save(i,j)=0.
           enddo
           rotzplot(i)=0.
           rotzplot_save(i)=0.
        enddo
C-- do the actual work
        open(38,file=fnaam)
        write(*,*)'Y and Z EPSILONS from file ',fnaam
 10     continue
        read(38,*,ERR=99,END=99)
     +          id,dy1,dy2,dy3,dz1,dz2,dz3,rot
C-- offer some protection against stupidities
        if (id.lt.0 .or. id.gt.100) then
           print *,'unknown id eps file ',id
*     goto 10 ! or stop
           stop
        endif
c        print *,'id dy dz rotplot'
c	print *,id,dy1,dy2,dy3,dz1,dz2,dz3,rot
        
        if (id.gt.0) then
        dyeps(id,1)=dy1
        dyeps(id,2)=dy2
        dyeps(id,3)=dy3
        dzeps(id,1)=dz1
        dzeps(id,2)=dz2
        dzeps(id,3)=dz3
        rotzplot(id)=rot
	endif

        
        goto 10
 99     continue
        close(38)



	end


	subroutine deltaz_final
        include 'Online.cmm'
        include 'getcor.cmm'
	common /kang/anw,ane,azw,aze,dzdy
	real dzdy(100)
        dy13=2.*26.011
	s=-1.
C--
C step 1 to 3
C--
C-- dznw
        dznw13=0.000 ! assumed
C-- angles 
        anw=(dzeps(1,3)-dzeps(1,1) - dznw13)/dy13
	ane=(dzeps(2,3)-dzeps(2,1) - dznw13)/dy13
	azw=s*(dzeps(5,3)-dzeps(5,1) - dznw13)/dy13
	aze=s*(dzeps(6,3)-dzeps(6,1) - dznw13)/dy13
C-- dzzw 1 to 3 	
	dzzw_p=dzeps(3,3)-dzeps(3,1) - azw*dy13
	dzze_p=dzeps(4,3)-dzeps(4,1) - aze*dy13
	dzzw_m=dzeps(7,3)-dzeps(7,1) - s*anw*dy13
	dzze_m=(dzeps(8,3)-dzeps(8,1) - s*ane*dy13)
	print *,'ASSUMED dz NW step 1 to 3: ',dznw13
	write(*,3002)anw,ane,azw,aze
 3002	format(' angle NW NE SW SE ',4F7.4)
	write(*,3001)dzzw_p,dzze_p,dzzw_m,dzze_m
 3001	format(' South: step 1 to 3 dz SW+ SE+ SW- SE- ',4F7.3)
C--
C step 1 to 2
C--
	il=1
	dy=26.011
        dznw_p=(dzeps(1,il+1)-dzeps(1,il) - anw*dy)
        dznw_m=(dzeps(5,il+1)-dzeps(5,il) - s*azw*dy)
        dzne_p=(dzeps(2,il+1)-dzeps(2,il) - ane*dy)
        dzne_m=(dzeps(6,il+1)-dzeps(6,il) - s*aze*dy)
	write(*,3003)dznw_p,dzne_p,dznw_m,dzne_m
 3003	format(' North: step 1 to 2 dz NW+ NE+ NW- NE- ',4F7.3)

        dzzw_p=dzeps(3,il+1)-dzeps(3,il) - azw*dy
        dzzw_m=dzeps(7,il+1)-dzeps(7,il) - s*anw*dy
        dzze_p=(dzeps(4,il+1)-dzeps(4,il) - aze*dy)
        dzze_m=(dzeps(8,il+1)-dzeps(8,il) - s*ane*dy)
	write(*,3004)dzzw_p,dzze_p,dzzw_m,dzze_m
 3004	format(' South: step 1 to 2 dz SW+ SE+ SW- SE- ',4F7.3)
C-- 
C step 2 to 3
C--
        il=2
        dy=26.011
        dznw_p=(dzeps(1,il+1)-dzeps(1,il) - anw*dy)
        dznw_m=(dzeps(5,il+1)-dzeps(5,il) - s*azw*dy)
        dzne_p=(dzeps(2,il+1)-dzeps(2,il) - ane*dy)
        dzne_m=(dzeps(6,il+1)-dzeps(6,il) - s*aze*dy)
        write(*,3005)dznw_p,dzne_p,dznw_m,dzne_m
 3005   format(' North: step 2 to 3 dz NW+ NE+ NW- NE- ',4F7.3)

        dzzw_p=dzeps(3,il+1)-dzeps(3,il) - azw*dy
        dzzw_m=dzeps(7,il+1)-dzeps(7,il) - s*anw*dy
        dzze_p=(dzeps(4,il+1)-dzeps(4,il) - aze*dy)
        dzze_m=(dzeps(8,il+1)-dzeps(8,il) - s*ane*dy)
        write(*,3006)dzzw_p,dzze_p,dzzw_m,dzze_m
 3006   format(' South: step 2 to 3 dz SW+ SE+ SW- SE- ',4F7.3)


       print *,'dy  NW NE SW SE with these angles'
	dzdy(1)=anw
	dzdy(2)=ane
	dzdy(3)=azw
	dzdy(4)=aze
	dzdy(5)=s*azw	
	dzdy(6)=s*aze	
	dzdy(7)=s*anw	
	dzdy(8)=s*ane	
        do id=1,8
           step12=(dyeps(id,2)-dyeps(id,1))-15.017*dzdy(id) !
           step23=(dyeps(id,3)-dyeps(id,2))+15.017*dzdy(id) !
           step13=dyeps(id,3)-dyeps(id,1)
C          print *,'id ',id,' alfa ',dzdy(id)
           write(*,5000)id,step12,step23,dzdy(id)*1000
 5000   format('id',1X,I2,' step 12 23  ',2F7.3,' used angle: ',F7.4)
        enddo
	end
