(* ::Package:: *)

(* Defining expressions for physical quantites *)

FractionUnfolded =
      Compile[{{midpoint, _Real}, {dHfold0, _Real}, {dCpfold, _Real}, {T,_Real}},
1/(1 + (E^((1/0.00831441)*( dHfold0*((-1/T) + (1/midpoint))
+ dCpfold *((midpoint/T) + Log[T/midpoint] - 1)))))];

CDobsUnfolding =
      Compile[{{cdF0, _Real}, {mF, _Real}, {cdU0, _Real}, {mU, _Real},
{midpoint, _Real}, {dHfold0, _Real}, {dCpfold, _Real}, {T, _Real}},
cdF0 - (mF*midpoint) + (mF*T) + (-cdF0 + cdU0 + (mF - mU)*(midpoint - T))/
(1 + ((E^(((midpoint - T)*(midpoint*dCpfold - dHfold0))/
(0.0083178*midpoint*T)))*((T/midpoint)^(dCpfold/0.0083178))))];

DeltaHobs =
      Compile[{{midpoint, _Real}, {dHfold0, _Real}, {dCpfold, _Real},
{dHbind0, _Real}, {dCpbind, _Real}, {Tref, _Real}, {T, _Real}},
((T - Tref)*dCpbind) + dHbind0 +
(((T*dCpfold) + dHfold0 - (midpoint*dCpfold))/
(1 + (E^(((((midpoint - T)*
(midpoint*dCpfold - dHfold0))/(midpoint*T)) + (dCpfold*Log[T/midpoint]))
/0.0083178))))];

DeltaGobs =
      Compile[{{midpoint, _Real}, {dHfold0, _Real}, {dCpfold, _Real},
{dHbind0, _Real}, {dCpbind, _Real}, {dSbind0, _Real}, {Tref, _Real}, {T,
_Real}}, T*dCpbind - Tref*dCpbind +
          dHbind0 - T*dSbind0 +
          0.0083178*T*
            Log[1 + (E^((dCpfold - midpoint*dCpfold/T -
dHfold0/midpoint + dHfold0/T - dCpfold*Log[T/midpoint])/0.0083178))] -
          T*dCpbind*Log[T/Tref]];

(* Fitting algorithm - LevenbergMarquardt alternating with parameter random modification *)

CoupledFoldingBindingFit[cdobs_, Ta_, errorcd_, dhobs_,  Tb_,
        errordhobs_, dgobs_, Tc_, errordgobs_, reftemp_] :=
      Module[{dp, dadp, dbdp, dcdp, \[Alpha], a, \[Beta], weightinga,
          weightingb, weightingc, ba, bb,
          bc, \[CapitalDelta]ia, \[CapitalDelta]ib, \[CapitalDelta]ic,
\[Chi]old, \[Chi], \[Lambda], parameterguess, p, \[CapitalDelta]p,
          loop, \[CapitalDelta]itesta, \[CapitalDelta]itestb,
\[CapitalDelta]itestc},

parameterguess = {Min[cdobs], 0.01, Max[cdobs], 0.01, Median[Ta],
dhobs[[Length[dhobs]]] - dhobs[[1]],
(dhobs[[Length[dhobs]]] - dhobs[[IntegerPart[Length[dhobs]/2]]])/
               (Tb[[Length[dhobs]]] - Tb[[IntegerPart[Length[dhobs]/2]]]),
dhobs[[1]],
(dhobs[[IntegerPart[Length[dhobs]/2]]] - dhobs[[1]])/
               (Tc[[IntegerPart[Length[dhobs]/2]]] - Tc[[1]]),
            0.01};
dp = 0.000000005;
dadp = Table[0.0, {10}, {Length[cdobs]}];
dbdp = Table[0.0, {10}, {Length[dhobs]}];
dcdp = Table[0.0, {10}, {Length[dgobs]}];
\[Alpha] = Table[0.0, {10}, {10}];
\[Beta] = Table[0.0, {10}];
If[((errorcd . errorcd) >
                0.0000001) && ((errordhobs . errordhobs) >
                0.0000001) && ((errordgobs . errordgobs) >
                0.0000001),
weightinga = 1.0/(errorcd*errorcd);
weightingb = 1.0/(errordhobs*errordhobs);
weightingc = 1.0/(errordgobs*errordgobs); ,
weightinga = Table[100.0, {Length[cdobs]}];
weightingb = Table[10.0, {Length[dhobs]}];
weightingc = Table[100.0, {Length[dgobs]}]; ];
p = parameterguess;
\[CapitalDelta]ia = cdobs - Table[
              CDobsUnfolding[p[[1]], p[[2]],
                p[[3]], p[[4]], p[[5]],
                p[[6]], p[[7]],
                Ta[[k]]], {k, 1,
                Length[cdobs]}];
\[CapitalDelta]ib = dhobs - Table[
              DeltaHobs[p[[5]], p[[6]],
                p[[7]], p[[8]], p[[9]],
                reftemp, Tb[[k]]], {k, 1,
                Length[dhobs]}];
\[CapitalDelta]ic = dgobs - Table[
              DeltaGobs[p[[5]], p[[6]],
                p[[7]], p[[8]], p[[9]],
                p[[10]], reftemp, Tc[[k]]], {k, 1,
                Length[dgobs]}];
bestp = parameterguess;
best\[Chi] =  (weightinga*\[CapitalDelta]ia) . \[CapitalDelta]ia +
              (weightingb*\[CapitalDelta]ib) . \[CapitalDelta]ib +
              (weightingc*\[CapitalDelta]ic) . \[CapitalDelta]ic;

Do[
  Do[
            p[[k]] *= RandomReal[{1.005, 0.995}]; , {k, 1,10}];

            \[CapitalDelta]ia =
            cdobs - Table[
            CDobsUnfolding[p[[1]], p[[2]],
                  p[[3]], p[[4]], p[[5]],
                  p[[6]], p[[7]],
                  Ta[[k]]], {k, 1,
                  Length[cdobs]}];

            \[CapitalDelta]ib =
            dhobs - Table[
                DeltaHobs[p[[5]], p[[6]],
                  p[[7]], p[[8]], p[[9]],
                  reftemp, Tb[[k]]], {k, 1,
                  Length[dhobs]}];
            \[CapitalDelta]ic =
            dgobs - Table[
                DeltaGobs[p[[5]], p[[6]],
                  p[[7]], p[[8]], p[[9]],
                  p[[10]], reftemp, Tc[[k]]], {k, 1,
                  Length[dgobs]}];

            \[Chi]old = (weightinga*\[CapitalDelta]ia) .  \[CapitalDelta]ia +
                        (weightingb*\[CapitalDelta]ib) .  \[CapitalDelta]ib +
                        (weightingc*\[CapitalDelta]ic) .  \[CapitalDelta]ic;
 Do[
            dadp[[1,k]] = (CDobsUnfolding[p[[1]] + dp,
                      p[[2]], p[[3]],
                      p[[4]], p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]] - dp,
                      p[[2]], p[[3]],
                      p[[4]], p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[2, k]] = (CDobsUnfolding[p[[1]],
                      p[[2]] + dp, p[[3]],
                      p[[4]], p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]],
                      p[[2]] - dp, p[[3]],
                      p[[4]], p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[3,k]] = (CDobsUnfolding[p[[1]],
                      p[[2]], p[[3]] + dp,
                      p[[4]], p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]], p[[2]],
                      p[[3]] - dp, p[[4]],
                      p[[5]], p[[6]],
                      p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[4,k]] = (CDobsUnfolding[p[[1]],
                      p[[2]], p[[3]],
                      p[[4]] + dp, p[[5]],
                      p[[6]], p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]], p[[2]],
                      p[[3]], p[[4]] - dp,
                      p[[5]], p[[6]],
                      p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[5,k]] = (CDobsUnfolding[p[[1]],
                      p[[2]], p[[3]],
                      p[[4]], p[[5]] + dp,
                      p[[6]], p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]], p[[2]],
                      p[[3]], p[[4]],
                      p[[5]] - dp, p[[6]],
                      p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[6,k]] = (CDobsUnfolding[p[[1]],
                      p[[2]], p[[3]],
                      p[[4]], p[[5]],
                      p[[6]] + dp, p[[7]] ,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]], p[[2]],
                      p[[3]], p[[4]],
                      p[[5]], p[[6]] - dp,
                      p[[7]] ,
                      Ta[[k]]])/(2*dp);
            dadp[[7,k]] = (CDobsUnfolding[p[[1]],
                      p[[2]], p[[3]],
                      p[[4]], p[[5]],
                      p[[6]], p[[7]]  + dp,
                      Ta[[k]]] -
                    CDobsUnfolding[p[[1]], p[[2]],
                      p[[3]], p[[4]],
                      p[[5]], p[[6]],
                      p[[7]] - dp ,
                      Ta[[k]]])/(2*dp);
            dadp[[8, k]] = 0.0;
            dadp[[9, k]] = 0.0;
            dadp[[10, k]] = 0.0; , {k, 1,Length[cdobs]}];

Do[
            dbdp[[1, k]] = 0.0;
            dbdp[[2, k]] = 0.0;
            dbdp[[3, k]] = 0.0;
            dbdp[[4, k]] = 0.0;
            dbdp[[5,k]] = (DeltaHobs[p[[5]] + dp,
                      p[[6]], p[[7]] ,
                      p[[8]], p[[9]], reftemp,
                      Tb[[k]]] -
                    DeltaHobs[p[[5]] - dp, p[[6]],
                      p[[7]] , p[[8]],
                      p[[9]], reftemp,
                      Tb[[k]]])/(2*dp);
            dbdp[[6,k]] = (DeltaHobs[p[[5]],
                      p[[6]] + dp, p[[7]] ,
                      p[[8]], p[[9]], reftemp,
                      Tb[[k]]] -
                    DeltaHobs[p[[5]], p[[6]] - dp,
                      p[[7]] , p[[8]],
                      p[[9]], reftemp,
                      Tb[[k]]])/(2*dp);
            dbdp[[7,k]] = (DeltaHobs[p[[5]],
                      p[[6]], p[[7]] + dp ,
                      p[[8]], p[[9]], reftemp,
                      Tb[[k]]] -
                    DeltaHobs[p[[5]], p[[6]],
                      p[[7]]  - dp, p[[8]],
                      p[[9]], reftemp,
                      Tb[[k]]])/(2*dp);
            dbdp[[8,k]] = (DeltaHobs[p[[5]],
                      p[[6]], p[[7]] ,
                      p[[8]] + dp, p[[9]], reftemp,
                      Tb[[k]]] -
                    DeltaHobs[p[[5]], p[[6]],
                      p[[7]] , p[[8]] - dp,
                      p[[9]], reftemp,
                      Tb[[k]]])/(2*dp);
            dbdp[[9,k]] = (DeltaHobs[p[[5]],
                      p[[6]], p[[7]] ,
                      p[[8]], p[[9]] + dp, reftemp,
                      Tb[[k]]] -
                    DeltaHobs[p[[5]], p[[6]],
                      p[[7]] , p[[8]],
                      p[[9]] - dp, reftemp,
                      Tb[[k]]])/(2*dp);
            dbdp[[10, k]] = 0.0; , {k, 1,
              Length[dhobs]}];
Do[
            dcdp[[1, k]] = 0.0;
            dcdp[[2, k]] = 0.0;
            dcdp[[3, k]] = 0.0;
            dcdp[[4, k]] = 0.0;
            dcdp[[5,k]] = (DeltaGobs[p[[5]] + dp,
                      p[[6]], p[[7]] ,
                      p[[8]], p[[9]],
                      p[[10]], reftemp, Tc[[k]]] -
                    DeltaGobs[p[[5]] - dp, p[[6]],
                      p[[7]] , p[[8]],
                      p[[9]], p[[10]], reftemp,
                      Tc[[k]]])/(2*dp);
            dcdp[[6,k]] = (DeltaGobs[p[[5]],
                      p[[6]] + dp, p[[7]] ,
                      p[[8]], p[[9]],
                      p[[10]], reftemp, Tc[[k]]] -
                    DeltaGobs[p[[5]], p[[6]] - dp,
                      p[[7]] , p[[8]],
                      p[[9]], p[[10]], reftemp,
                      Tc[[k]]])/(2*dp);
            dcdp[[7,k]] = (DeltaGobs[p[[5]],
                      p[[6]], p[[7]] + dp ,
                      p[[8]], p[[9]],
                      p[[10]], reftemp, Tc[[k]]] -
                    DeltaGobs[p[[5]], p[[6]],
                      p[[7]]  - dp, p[[8]],
                      p[[9]], p[[10]], reftemp,
                      Tc[[k]]])/(2*dp);
            dcdp[[8,k]] = (DeltaGobs[p[[5]],
                      p[[6]], p[[7]] ,
                      p[[8]] + dp, p[[9]],
                      p[[10]], reftemp, Tc[[k]]] -
                    DeltaGobs[p[[5]], p[[6]],
                      p[[7]] , p[[8]] - dp,
                      p[[9]], p[[10]], reftemp,
                      Tc[[k]]])/(2*dp);
            dcdp[[9,k]] = (DeltaGobs[p[[5]],
                      p[[6]], p[[7]] ,
                      p[[8]], p[[9]] + dp,
                      p[[10]], reftemp, Tc[[k]]] -
                    DeltaGobs[p[[5]], p[[6]],
                      p[[7]] , p[[8]],
                      p[[9]] - dp, p[[10]], reftemp,
                      Tc[[k]]])/(2*dp);
            dcdp[[10,k]] = (DeltaGobs[p[[5]],
                      p[[6]], p[[7]] ,
                      p[[8]], p[[9]],
                      p[[10]] + dp, reftemp,
                      Tc[[k]]] -
                    DeltaGobs[p[[5]], p[[6]],
                      p[[7]] , p[[8]],
                      p[[9]], p[[10]] - dp, reftemp,
                      Tc[[k]]])/(2*dp); , {k, 1,Length[dgobs]}];


              \[Alpha] = \[Alpha]*0.0;
              Do[Do[Do[\[Alpha][[i, j]] +=
                  weightinga[[k]]*dadp[[i, k]]*dadp[[j, k]]; , {k, 1,Length[cdobs]}];
                  Do[\[Alpha][[i,j]] +=
                  weightingb[[k]]*dbdp[[i, k]]*dbdp[[j, k]]; , {k, 1,Length[dhobs]}];
                  Do[\[Alpha][[i,j]] += weightingc[[k]]*dcdp[[i, k]]* dcdp[[j, k]]; ,
                  {k, 1,Length[dgobs]}]; , {i, 1,10}]; , {j, 1,10}];

              loop =0;
              \[Lambda] = 0.1;
              \[CapitalDelta]p =
            Table[1.0, {10}]; While[(Apply[Plus,
                    Abs[\[CapitalDelta]p]] > 0.00001) &&
              loop < 101, loop += 1;
ba = dadp . (weightinga* \[CapitalDelta]ia);
bb = dbdp . (weightingb* \[CapitalDelta]ib);
bc = dcdp . (weightingc* \[CapitalDelta]ic);
\[Beta] = ba + bb + bc;
a = \[Alpha];
Do[
              a[[k,k]] = (1 + \[Lambda])*\[Alpha][[k,k]], {k, 1,10}];
              \[CapitalDelta]p = Inverse[a] . \[Beta];
              \[CapitalDelta]itesta = cdobs - Table[
                  CDobsUnfolding[
                    p[[1]] + \[CapitalDelta]p[[1]],
                    p[[2]] + \[CapitalDelta]p[[2]],
                    p[[3]] + \[CapitalDelta]p[[3]],
                    p[[4]] + \[CapitalDelta]p[[4]],
                    p[[5]] + \[CapitalDelta]p[[5]],
                    p[[6]] + \[CapitalDelta]p[[6]],
                    p[[7]] + \[CapitalDelta]p[[7]],
                    Ta[[k]]], {k, 1,
                    Length[cdobs]}];

              \[CapitalDelta]itestb = dhobs - Table[DeltaHobs[
                    p[[5]] + \[CapitalDelta]p[[5]],
                    p[[6]] + \[CapitalDelta]p[[6]],
                    p[[7]] + \[CapitalDelta]p[[7]],
                    p[[8]] + \[CapitalDelta]p[[8]],
                    p[[9]] + \[CapitalDelta]p[[9]],
                    reftemp, Tb[[k]]], {k, 1, Length[dhobs]}];

              \[CapitalDelta]itestc = dgobs - Table[DeltaGobs[
                    p[[5]] + \[CapitalDelta]p[[5]],
                    p[[6]] + \[CapitalDelta]p[[6]],
                    p[[7]] + \[CapitalDelta]p[[7]],
                    p[[8]] + \[CapitalDelta]p[[8]],
                    p[[9]] + \[CapitalDelta]p[[9]],
                    p[[10]] + \[CapitalDelta]p[[10]],
                    reftemp, Tc[[k]]], {k, 1, Length[dgobs]}];

                    \[Chi] =
(weightinga*\[CapitalDelta]itesta) .  \[CapitalDelta]itesta +
(weightingb*\[CapitalDelta]itestb) .  \[CapitalDelta]itestb +
(weightingc*\[CapitalDelta]itestc) .  \[CapitalDelta]itestc;
            If[0.999999999999 > (\[Chi]/\[Chi]old), \[Lambda] *=
                0.1; \[Chi]old = \[Chi]; \[CapitalDelta]ia =
\[CapitalDelta]itesta; \[CapitalDelta]ib = \[CapitalDelta]itestc;
\[CapitalDelta]ic = \[CapitalDelta]itestc; p += \[CapitalDelta]p; ];
            If[1.0 < (\[Chi]/\[Chi]old), \[Lambda] *=
                10.0; ]; ];
If[\[Chi]old < best\[Chi], best\[Chi] = \[Chi]old;
            bestp = p; ]; , {l, 1,200}];

Return[{bestp, best\[Chi]}];];

