diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp
index 7eb8507..2f61884 100644
--- a/physics/2pt-fit.cpp
+++ b/physics/2pt-fit.cpp
@@ -16,7 +16,7 @@ int main(int argc, char *argv[])
     // parse arguments /////////////////////////////////////////////////////////
     OptParser            opt;
     bool                 parsed, doPlot, doHeatmap, doCorr, fold;
-    string               corrFileName, model, outFileName, outFmt;
+    string               corrFileName, model, outFileName, outFmt, save;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
     Minimizer::Verbosity verbosity;
@@ -30,7 +30,7 @@ int main(int argc, char *argv[])
     opt.addOption("s", "shift"    , OptParser::OptType::value  , true,
                   "time variable shift", "0");
     opt.addOption("m", "model"    , OptParser::OptType::value  , true,
-                  "fit model (exp|exp2|exp3|cosh|cosh2|cosh3|explin|const|<interpreter code>)", "cosh");
+                  "fit model (exp|exp2|exp3|sinh|cosh|cosh2|cosh3|cotanh|explin|const|<interpreter code>)", "cosh");
     opt.addOption("" , "nPar"     , OptParser::OptType::value  , true,
                   "number of model parameters for custom models "
                   "(-1 if irrelevant)", "-1");
@@ -48,6 +48,8 @@ int main(int argc, char *argv[])
                   "show the fit plot");
     opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
                   "show the fit correlation heatmap");
+    opt.addOption("s", "save", OptParser::OptType::value, true,
+                    "saves the source and .pdf", "");
     opt.addOption("", "help"      , OptParser::OptType::trigger, true,
                   "show this help message and exit");
     parsed = opt.parse(argc, argv);
@@ -71,6 +73,7 @@ int main(int argc, char *argv[])
     fold         = opt.gotOption("fold");
     doPlot       = opt.gotOption("p");
     doHeatmap    = opt.gotOption("h");
+    save         = opt.optionValue("s");
     switch (opt.optionValue<unsigned int>("v"))
     {
         case 0:
@@ -117,7 +120,7 @@ int main(int argc, char *argv[])
     
     // make model //////////////////////////////////////////////////////////////
     DoubleModel mod;
-    bool        coshModel = false, linearModel = false, constModel = false;
+    bool        sinhModel = false, coshModel = false, cotanhModel = false, linearModel = false, constModel = false;
     
     if ((model == "exp") or (model == "exp1"))
     {
@@ -144,6 +147,15 @@ int main(int argc, char *argv[])
                                    + p[5]*exp(-p[4]*x[0]);
                         }, 1, nPar);
     }
+    else if (model == "sinh")
+    {
+        sinhModel = true;
+        nPar      = 2;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[1]*(exp(-p[0]*x[0])-exp(-p[0]*(nt-x[0])));
+                        }, 1, nPar);
+    }
     else if ((model == "cosh") or (model == "cosh1"))
     {
         coshModel = true;
@@ -174,6 +186,15 @@ int main(int argc, char *argv[])
                                  + p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0])));
                         }, 1, nPar);
     }
+    else if (model == "cotanh")
+    {
+        cotanhModel = true;
+        nPar        = 2;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[0] + p[1]/tanh(0.1892*(nt/2-x[0])); // 0.1892 is the meson eff mass in lattice unit from 2pt fit
+                        }, 1, nPar);
+    }
     else if (model == "explin")
     {
         linearModel = true;
@@ -247,6 +268,11 @@ int main(int argc, char *argv[])
         init(0) = data.y(nt/4, 0)[central];
 
     }
+    else if(cotanhModel)
+    {
+        init(0) = 0.5*(data.y(0, 0)[central]+ data.y(nt-1, 0)[central]);
+        init(1) = 0.5*(data.y(0, 0)[central]- data.y(nt-1, 0)[central]);
+     }
     else
     {
         init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]);
@@ -269,9 +295,14 @@ int main(int argc, char *argv[])
         {
             globMin.setLowLimit(p, -10*fabs(init(0)));
             locMin.setLowLimit(p, -10*fabs(init(0)));
-            // cout << "Suppressing low limits" << endl;
             globMin.setHighLimit(p, 10*fabs(init(0)));
         }
+        else if(cotanhModel)
+        {
+            globMin.setLowLimit(p,-0.5);
+            locMin.setLowLimit(p, -0.3);
+            globMin.setHighLimit(p,1);
+        }
         else
         {
             globMin.setLowLimit(p, 0.);
@@ -328,6 +359,11 @@ int main(int argc, char *argv[])
         p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1);
         p << Color("rgb 'red'")  << PlotData(data.getData());
         p.display();
+        if(save != "")
+        {
+            cout << "Saving plot and source code to " << save << endl;
+            p.save(save);
+        }
         // effective mass plot //////////////////////////////////////////////////////
         if (!constModel)
         {