#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

//
// Functions to calculate R2eff as a function of
// exchange parameters and number of CPMG cycles.
//
// Baldwin et al. JMR 244, p114-124 (2014)
//
// Note: solution is exact for 2x2 Bloch McConnell
// equation. This treatment neglects:
// 1) Off resonance effects
// 2) Scalar coupling.
// 3) Different relaxation of different terms
// Each of which can lead to significant errors.
// For getting optimally accurate parameters from
// fitting data, it is advisable to use more 
// complete software such as CATIA.
//
// Compile with a c compiler, eg gcc.
//
//
// Copyright (c) A. Baldwin University of Oxford
// 15th March 2014


double cpmg_baldwin(double kex,double pb,double dw,double ncyc,double Trel,double R2g,double R2e)
{
  double pa=(1-pb);
  double keg=kex*(1-pb);
  double kge=kex*pb;
  double deltaR2=R2e-R2g;
  double nu_cpmg=ncyc/Trel;
  double tcp=Trel/(4.0*ncyc); // #time for one free precession element
  
  //#########################################################################
  //#get the real and imaginary components of the exchange induced shift
  double g1=2*dw*(deltaR2+keg-kge);                   //#same as carver richards zeta
  double g2=pow(deltaR2+keg-kge,2.)+4*keg*kge-dw*dw;   //#same as carver richards psi
  double g3=cos(0.5*atan2(g1,g2))*pow(g1*g1+g2*g2,0.25);   //#trig faster than square roots
  double g4=sin(0.5*atan2(g1,g2))*pow(g1*g1+g2*g2,0.25);   //#trig faster than square roots
  //#########################################################################
  //#time independent factors
  double complex N=(kge+g3-kge+I*g4);          //  #N=oG+oE
  double NNc=(g3*g3+g4*g4);
  double f0=(dw*dw+g3*g3)/(NNc);          //    #f0
  double f2=(dw*dw-g4*g4)/(NNc);           //   #f2
  //#t1=(-dw+g4)*(complex(-dw,-g3))/(NNc) #t1
  double complex t2=(dw+g4)*((dw-g3*I))/(NNc); //#t2
  double complex t1pt2=(2*dw*dw-g1*I)/(NNc);   //  #t1+t2
  double complex oGt2=((deltaR2+keg-kge-g3)+I*(dw-g4))*t2; // #-2*oG*t2
  
  double Rpre=(R2g+R2e+kex)/2.0;  // #-1/Trel*log(LpreDyn)
  
  //#do calc in numpy
  double E0= 2.0*tcp*g3; // #derived from relaxation       #E0=-2.0*tcp*(f00R-f11R)
  double E2= 2.0*tcp*g4; // #derived from chemical shifts  #E2=complex(0,-2.0*tcp*(f00I-f11I))
  double complex E1=((g3-g4*I))*tcp;  //  #mixed term (complex) (E0-iE2)/2
  double ex0b=(f0*cosh(E0)-f2*cos(E2));            //   #real
  double complex ex0c=(f0*csinh(E0)-f2*csin(E2)*I);   // #complex
  double complex ex1c=csinh(E1);                                  // #complex
  double complex v3=csqrt(ex0b*ex0b-1); //  #exact result for v2v3
  double complex y=cpow((ex0b-v3)/(ex0b+v3),ncyc);
  double complex v2pPdN=(( (deltaR2+kex+I*dw) )*ex0c+(-oGt2-kge*t1pt2)*2*ex1c);     //   #off diagonal common factor. sinh fuctions
  
  double complex Tog=(((1+y)/2+(1-y)/(2*v3)*(v2pPdN)/N));     
  double R2eff=Rpre-ncyc/(Trel)*acosh(creal(ex0b))-1/Trel*log((creal(Tog))); // #estimate R2eff
  return R2eff;
}
 

int main()
{

  double dfrq=200.;   //spectrometer frequency of nucleci (MHz) 
  double R2g=10.;     //relaxation rate of ground
  double R2e=10.;     //relaxation rate of excited
  double Trelax=0.02; //total time of CPMG block
  double kex=1000.;   //exchange rate (s-1)
  double pb=0.1;      //population of minor state
  double dw=1000;     //omegaE-omegaG (ppm)
  
  double pi=acos(-1);

  int ncycVals=19;     //number of ncyc points
  double ncyc_array[19]= {2,4,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40}; //#number of cpmg cycles
  
  for(int i=0;i<ncycVals;i++)
    {
      double ncyc=ncyc_array[i];
      double R2eff=cpmg_baldwin(kex,pb,dw*2*pi*dfrq,ncyc_array[i],Trelax,R2g,R2e);
    }
  
  return 0;
}