5 gemNULL, gemNONE, gemDD, gemAD, gemAA, gemA4, gemNR
7 static const char *gemType[] = {NULL, "none", "dd", "ad", "aa", "a4", NULL};
9 /* The first few sections of this file contain functions that were adopted,
10 * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
11 * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
12 * This is also the case with the function eq10v2() in geminate.c.
14 * The parts menetioned in the previous paragraph were contributed under a BSD license.
17 /* This first part is derived from complex.c which I recieved from Omer Markowitch.
20 * ------------- from complex.c ------------- */
23 /* definition of PI */
25 #define PI (acos(-1.0))
28 /* definition of the type complex */
35 /* ------------ end of complex.c ------------ */
37 /* This next part was derived from cerror.c and rerror.c,
38 * also received from Omer Markovitch.
39 * ------------- from [cr]error.c ------------- */
42 #define sPI (sqrt(PI))
45 /* ------------ end of [cr]error.c ------------ */
47 /* ///////////////// REVERSIBLE GEMINATE RECOMBINATION ///////////////////
48 * Here follow routines and structs for reversible geminate recombination.
60 /* Used in the rewritten version of Omer's gem. recomb. analysis */
61 double ka, kd; /* -- || -- results[] */
62 double sigma; /* -- || -- sigma */
63 double D; /* The diffusion coefficient */
64 double kD; /* Replaces kD in analyse_corr_gem3d() */
66 /* The following members are for calcsquare() and takeAwayBallistic() */
67 double tDelta; /* Time between frames */
68 /* double logAfterTime; /\* Time after which we do the lsq calculations on a logarithmic timescale. *\/ */
69 int nFitPoints; /* Number of points to take from the ACF for fitting */
70 double begFit; /* Fit from this time (ps) */
71 double endFit; /* Fit up to this time (ps) */
72 /* double logDelta; */
74 /* To get an equal number of points in the lin and log regime,
75 * we'll use logDelta to calculate where to sample the ACF.
76 * if i and j are indices in the linear and log regime, then:
77 * j = Cexp(A(i+nLin)),
78 * where C = (nLin**2 / len) and A = log(len/nLin) / nLin.
79 * This expands to j = (nLin**2 / len) * exp((i+nLin) * log(len/nLin) / nLin).
80 * We'll save part of our brains and some computing time if we pre-calculate
81 * 1) logDelta = log(len/nLin) / nLin
82 * 2) logPF = nLin**2 / len
83 * and get j = logPF * exp((i+nLin) * logDelta). */
85 /* I will redo this for a fit done entirely in log-log.
87 * nFitPoints' = nFitPoints-1
92 * (j'=len <=> i=nFitPoints')
93 * => A=log(len)/nFitPoints'
94 * => j = exp(i*log(len)/(nFitPoints-1)) -1
96 /* #define GETLOGINDEX(i,params) (params)->logPF * exp(((i)+(params)->nLin) * (params)->logDelta)
99 int nLin; /* Number of timepoints in the linear regime */
100 int len; /* Length of time and ct arrays */
101 int nExpFit; /* Number of exponentials to fit */
102 real ballistic; /* Time before which the ballistic term should be fitted */
103 gmx_bool bDt; /* TRUE => use time derivative at time 0
104 * to find fastest component.
105 * FALSE => use coefficient in exponenetial
106 * to find fastest component. */
111 size_t n; /* Number of data points (lin-log) */
112 double *y; /* The datapoints */
113 double *ctTheory; /* The theoretical ct which will be built by gemFunc_f. */
118 double tDelta; /* time difference between subsequent datapoints */
119 size_t nData; /* real size of the data */
121 double *doubleLogTime;
125 extern void takeAwayBallistic(double *ct, double *t,
127 int nexp, gmx_bool bDerivative);
130 extern t_gemParams *init_gemParams(const double sigma, const double D,
131 const real *t, const int len, const int nFitPoints,
132 const real begFit, const real endFit,
133 const real ballistic, const int nBalExp, const gmx_bool bDt);
135 /* Fit to geminate recombination model.
136 Returns root mean square error of fit. */
137 extern real fitGemRecomb(double *ct, double *time, double **ctFit,
138 const int nData, t_gemParams *params);
140 extern void dumpN(const real *e, const int nn, char *fn);
142 /* Fix NaN that might appear in the theoretical acf. */
143 extern void fixGemACF(double *ct, int len);