program spincoordinate
    use precision_def
    use sim_utils
    use compile_interface
    use muon_mod
  implicit none
!   type deriv_struct
!    real(rp) coef(4), coef_err(4), chisq
!    integer norder
!   end type
   type (deriv_struct) dspin_dx, dspin_dxp, dx_dp, dxp_dp, dp_dt
   character*100 filename
   real(rp) coef(4), coef_err(4), chisq
   real(rp) average_energy, stdev_energy
   real(rp) dgamma_dt, gamma_1, gamma_2
   real(rp) tau/2.197/,gamma_0/29.3/, t_rev/0.1492/  
   real(rp) Time,t  !microseconds
   real(rp) delta_1, delta_2, delta_omega_0, delta_omega_1, delta_omega_2, delta_omega_avg, delta_omega
   real(rp) omega_a
   real(rp) do0_err, do1_err, do2_err
   real(rp) delta_omega_avg_err
   integer lun
   integer ix, i
   character*200 string, directory
  
string='pwd > dir.dat'
call execute_command_line(string)
lun=lunget()
open(unit=lun,file='dir.dat')
read(lun,'(a)')string
i=1
do while(string(i:i)/='')
 if(string(i:i) == '/')ix=i
 i=i+1
end do
directory=string(ix+1:)
print *,'directory= ', directory

  filename = '2D_hist_fit_coef_err_7.dat'  
  call read_bin(filename, dspin_dx%coef, dspin_dx%coef_err, dspin_dx%chisq)
  filename = '2D_hist_fit_coef_err_8.dat'  
  call read_bin(filename, dspin_dxp%coef, dspin_dxp%coef_err, dspin_dxp%chisq)
  filename = '2D_hist_fit_coef_err_4.dat'  
  call read_bin(filename, dx_dp%coef, dx_dp%coef_err, dx_dp%chisq)
  filename = '2D_hist_fit_coef_err_5.dat'  
  call read_bin(filename, dxp_dp%coef, dxp_dp%coef_err, dxp_dp%chisq)
  filename='2D_hist_4.dat'
  call read_2D_hist(filename, average_energy, stdev_energy)
  dgamma_dt = stdev_energy**2/(1+average_energy)**2/tau
  print '(/,a,es12.4,a4)',' dgamma_dt=',dgamma_dt,' /us'
  delta_1 = average_energy
  delta_2 = dgamma_dt/gamma_0
!  omega_a = amu*gamma * omega_c = omega_c*gamma/(gamma**2-1)
   omega_a = twopi/t_rev * gamma_0/(gamma_0**2-1)
   print '(a,es12.4,a)','omega_a = ',omega_a,'rad/us'

   delta_omega_0 = dspin_dx%coef(2)*(dx_dp%coef(2)+2*dx_dp%coef(3)*delta_1 + 3*dx_dp%coef(4)*delta_1**2) + &
                   dspin_dxp%coef(2)*(dxp_dp%coef(2)+2*dxp_dp%coef(3)*delta_1 + 3*dxp_dp%coef(4)*delta_1**2)
   
   delta_omega_1 = dspin_dx%coef(2)*(2*dx_dp%coef(3)*delta_2 + 6*dx_dp%coef(4)*delta_1*delta_2) + &
                   dspin_dxp%coef(2)*(2*dxp_dp%coef(3)*delta_2 + 6*dxp_dp%coef(4)*delta_1*delta_2)

   delta_omega_2 = dspin_dx%coef(2)*3*dx_dp%coef(4)*delta_2**2  + dspin_dxp%coef(2)*3*dxp_dp%coef(4)*delta_2**2

   call get_errors(dspin_dx, dspin_dxp,dx_dp, dxp_dp,delta_1, delta_2, do0_err, do1_err, do2_err)

   t=0
   delta_omega = (delta_omega_0 + delta_omega_1 * t + delta_omega_2 * t**2) * delta_2

   Time = 600.
   delta_omega_avg = (delta_omega_0 + delta_omega_1 * Time/2. + delta_omega_2 * Time**2/3) * delta_2
   delta_omega_avg_err = sqrt((do0_err)**2 + (do1_err * Time/2.)**2 + (do2_err * Time**2/3)**2) * delta_2


   print '(2(a,es12.4))',' delta_1=',delta_1,' delta_2=',delta_2
   print '(3(a,es12.4,a3,es12.4))',' delta_omega_0 =',delta_omega_0 ,'+/-',do0_err,' delta_omega_1 =',delta_omega_1,'+/-',do1_err ,' delta_omega_2 =',delta_omega_2,'+/-',do2_err 
   print '(2(a,es12.4))',' delta_omega_avg=',delta_omega_avg,' delta_omega_avg/omega_a=',delta_omega_avg/omega_a
   print '(2(a,es12.4))',' delta_omega_avg_err=',delta_omega_avg_err,' delta_omega_avg_err/omega_a=',delta_omega_avg_err/omega_a
   write(6,'(2(a,es12.4))')'delta_omega_0=',delta_omega_0,'; do0_err=',do0_err
   write(6,'(2(a,es12.4))')'delta_omega_1=',delta_omega_1,'; do1_err=',do1_err
   write(6,'(2(a,es12.4))')'delta_omega_2=',delta_omega_2,'; do2_err=',do2_err
   write(6,'(2(a,es12.4))')'delta_2=',delta_2,'; delta_1=',delta_1
   write(6,'(a,es12.4)')'omega_a=',omega_a

   lun=lunget()
   open(unit=lun,file='domega_time_dependence.dat')
   write(lun,'(2(a,es12.4))')'delta_omega_0=',delta_omega_0,'; do0_err=',do0_err
   write(lun,'(2(a,es12.4))')'delta_omega_1=',delta_omega_1,'; do1_err=',do1_err
   write(lun,'(2(a,es12.4))')'delta_omega_2=',delta_omega_2,'; do2_err=',do2_err
   write(lun, '(3(a,es12.4,a11,es12.4))')' delta_omega_0 =',delta_omega_0 ,'; do0_err=',do0_err,'; delta_omega_1 =',delta_omega_1,'; do1_err= ',do1_err ,'; delta_omega_2 =',delta_omega_2,'; do2_err= ',do2_err 
   write(lun, '(2(a,es12.4))')' delta_omega_avg=',delta_omega_avg,'; delta_omega_avg_omega_a=',delta_omega_avg/omega_a
   write(lun, '(2(a,es12.4))')' delta_omega_avg_err=',delta_omega_avg_err,'; delta_omega_avg_err_omega_a=',delta_omega_avg_err/omega_a

   write(lun,'(2(a,es12.4))')'delta_2=',delta_2,'; delta_1=',delta_1
   write(lun,'(a,es12.4)')'omega_a=',omega_a
   close(unit=lun)
   print '(a)',' write domega_time_dependence.dat'
   lun=lunget()
   open(unit=lun, file='../spin_summary.dat',access='append')
   write(lun, '(a12,2(a20),3(a30),4a20)')'# directory','delta_omega_avg','delta_omega_avg_err','delta_omega_avg/omega_a','delta_omega_avg_err/omega','dphi/dx','dphi/dx err','dphi/dxp','dphi/dxp err'
   write(lun, '(a12,2es20.4,3es30.4,4es20.4)')trim(directory),delta_omega_avg, delta_omega_avg_err,delta_omega_avg/omega_a,delta_omega_avg_err/omega_a,dspin_dx%coef(2),dspin_dx%coef_err(2),dspin_dxp%coef(2),dspin_dxp%coef_err(2)



  end program spincoordinate

  subroutine read_bin(filename, coef, coef_err, chisq)
   use sim_utils
   implicit none
   integer lun, reason
   integer number
   integer ix, ix_equal, ix_semi,i,j
   real(rp) coef(4), coef_err(4), chisq
   character*1 cn(10)/'0','1','2','3','4','5','6','7','8','9'/
   character*100, string, filename, string1, string2, string3
   character*10 word
   lun=lunget()


   open(unit=lun,file=filename)
    
    read(lun, '(a)')string
      if(index(string,'number')/=0)read(string(index(string,'=')+1:),*)number
      print '(a,i10)',' number = ',number
      read(lun,'(a)')string
!      print '(a)',string
      i=0
      do while (.true.)
        i=i+1
        word = trim('coef'//cn(i))
        ix = index(string,trim(word))
        if(ix == 0)exit
        string2=string(ix+5:)
        ix_equal=index(string2,'=')
        string3=string2(ix_equal+1:)
        read(string3, * )coef(i)
      end do
      print '(a,4es12.4)','coef(j) =',(coef(j),j=1,i-1)

      read(lun,'(a)')string
!      print '(a)',string
      i=0
      do while (.true.)
        i=i+1
        word = trim('coef'//cn(i))//'_err'
        ix = index(string,trim(word))
        if(ix == 0)exit
        string2=string(ix+5:)
        ix_equal=index(string2,'=')
        string3=string2(ix_equal+1:)
        read(string3, * )coef_err(i)
!        print '(a,3x,i10,3x,a,es12.4)','i=',i,'coef_err(i) =',coef_err(i)
      end do
      print '(a,4es12.4)','coef_err(j) =',(coef_err(j),j=1,i-1)    
      read(lun,'(a)')string
      if(index(string,'chisq')/=0)read(string(index(string,'=')+1:),*)chisq
      print '(a,es12.4,/)','chisq=',chisq
    close(lun)    
  return
end subroutine
        
subroutine read_2D_hist(filename, average_energy, stdev_energy)
   use sim_utils
   implicit none
   real(rp) average_energy, stdev_energy
   character*100, string, filename
   integer lun
   integer ix, ix1

   lun=lunget()
   open(unit=lun, file=filename)
   do while(.true.)
    read(lun,'(a)')string
    ix= index(string, 'average_energy')
    if(ix == 0)cycle
      ix=index(string,'=')
      read(string(ix+1:),*)average_energy
      ix1 = index(string(ix+2:),'=')
      read(string(ix+2+ix1+1:),*)stdev_energy
      exit
   end do
   print '(a,es12.4,a,es12.4)','average_energy=', average_energy,' stdev_energy=',stdev_energy
   return
  end subroutine read_2D_hist
   
  subroutine get_errors(dspin_dx, dspin_dxp,dx_dp, dxp_dp,delta1, delta2, do0_err, do1_err, do2_err)
    use precision_def
    use sim_utils
    use muon_mod
!    use compile_interface
  implicit none
   type (deriv_struct) dspin_dx, dspin_dxp, dx_dp, dxp_dp
   real(rp) delta1, delta2
   real(rp) A,B,C,D,E,F,G,H,J,K
   real(rp) do0_err, do1_err, do2_err
   real(rp) do0_dA,do0_dB,do0_dC,do0_dD,do0_dE,do0_dF,do0_dG,do0_dH,do0_dJ
   real(rp) do1_dA,do1_dC,do1_dD,do1_dE,do1_dH,do1_dJ, do1_dK
   real(rp) do2_dA,do2_dE,do2_dK,do2_dF,do2_dJ
   real(rp) A_err,B_err,C_err,D_err,E_err,F_err,G_err,H_err,J_err,K_err
   real(rp) delta1_err, delta2_err
   
  delta1_err=0.
  delta2_err=0.
  A=  dspin_dx%coef(2); B=dx_dp%coef(2); C=dx_dp%coef(3); D=delta1;E= dx_dp%coef(4);F= dspin_dxp%coef(2);G=dxp_dp%coef(2);H=dxp_dp%coef(3);J=dxp_dp%coef(4)
  A_err=dspin_dx%coef_err(2); B_err=dx_dp%coef_err(2); C_err=dx_dp%coef_err(3); D_err=delta1_err
  E_err= dx_dp%coef_err(4);F_err= dspin_dxp%coef_err(2);G_err=dxp_dp%coef_err(2);H_err=dxp_dp%coef_err(3);J_err=dxp_dp%coef_err(4)
  K=delta2
  K_err=delta2_err

!   delta_omega_0 = A*(B+2*C*D + 3*E*D**2) + F*(G+2*H*D + 3*J*D**2)
   
!   delta_omega_1 = A*(2*C*K + 6*E*D*K) + F*(2*H*K + 6*J*D*K)

!   delta_omega_2 = A*3*E*K**2  + F*3*J*K**2

   do0_dA = (B+2*C*D + 3*E*D**2)
   do0_dB = A
   do0_dC = 2*A*D
   do0_dD = A*(2*C+6*E*D)+ F*(2*H + 6*J*D)
   do0_dE = 3*A*D**2
   do0_dF = (G+2*H*D + 3*J*D**2)
   do0_dG= F
   do0_dH=2*F*D
   do0_dJ= 3*F*D**2
   do0_err =sqrt( (do0_dA*A_err)**2+(do0_dB*B_err)**2+(do0_dC*C_err)**2+(do0_dD*D_err)**2+(do0_dE*E_err)**2+(do0_dF*F_err)**2+(do0_dG*G_err)**2+(do0_dH*H_err)**2+(do0_dJ*J_err)**2)

   do1_dA =(2*C*K + 6*E*D*K)
   do1_dC=2*A*K
   do1_dD=6*A*E*K + 6*F*J*K
   do1_dE=6*A*D*K
   do1_dK = A*(2*C+6*E*D)+F*(2*H+6*J*D)
   do1_dH= 2*F*K
   do1_dJ=6*F*D*K
   do1_err=sqrt((do1_dA*A_err)**2+(do1_dC*C_err)**2+(do1_dD*D_err)**2+(do1_dE*E_err)**2+(do1_dK*K_err)**2+(do1_dH*H_err)**2+(do1_dJ*J_err)**2)

   do2_dA=3*E*K**2
   do2_dE=3*A*K**2
   do2_dK=2*(A*3*E*K  + F*3*J*K)
   do2_dF = 3*J*K**2
   do2_dJ=3*F*K**2
   do2_err=sqrt((do2_dA*A_err)**2+(do2_dE*E_err)**2+(do2_dK*K_err)**2+(do2_dF*F_err)**2+(do2_dJ*J_err)**2)

   return
   end subroutine get_errors