Changeset 20936 for trunk/Ohana/src/opihi/cmd.data/vgauss.c

Timestamp:

Dec 7, 2008, 3:31:01 PM (18 years ago)

Author:

eugene

Message:

big update from eam_branch_20081124 with updates to Opihi math

File:

• trunk/Ohana/src/opihi/cmd.data/vgauss.c (modified) (7 diffs)

trunk/Ohana/src/opihi/cmd.data/vgauss.c

@@ -2 +2 @@
 
 /* local private functions */
-float fgaussOD (float, float *, int, float *);
+opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
 
 # define GET_VAR(V,A) \
@@ -16 +16 @@
 
   int i, N, Npts, Npar, Quiet;
-  float par[4], *v1, *v2, *dy, **covar;
-  float chisq, ochisq, dchisq;
+  opihi_flt par[4], *v1, *v2, *dy, **covar;
+  opihi_flt chisq, ochisq, dchisq;
   Vector *xvec, *yvec, *svec, *ovec;
   char *c, name[16];
@@ -42 +42 @@
   if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
+  REQUIRE_VECTOR_FLT (xvec, FALSE); 
+  REQUIRE_VECTOR_FLT (yvec, FALSE); 
+  REQUIRE_VECTOR_FLT (svec, FALSE); 
+
   Npts = xvec[0].Nelements;
-  ALLOCATE (dy, float, Npts);
-  REALLOCATE (ovec[0].elements, float, Npts);
+  ResetVector (ovec, OPIHI_FLT, Npts);
+
+  ALLOCATE (dy, opihi_flt, Npts);
 
   GET_VAR (par[0], "C0");
@@ -52 +57 @@
   Npar = 4;
 
-  v1 = svec[0].elements;
+  // mrqmin takes the inverse variance (do not generate NANs)
+  v1 = svec[0].elements.Flt;
   v2 = dy;
   for (i = 0; i < Npts; i++, v1++, v2++) {
@@ -58 +64 @@
   } 
   
-  ochisq = mrqinit (xvec[0].elements, yvec[0].elements, dy, Npts, par, Npar, fgaussOD, !Quiet);
+  ochisq = mrqinit (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, fgaussOD, !Quiet);
   dchisq = ochisq + 2*Npts;
 
   for (i = 0; (i < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); i++) {
-    chisq = mrqmin (xvec[0].elements, yvec[0].elements, dy, Npts, par, Npar, fgaussOD, !Quiet);
+    chisq = mrqmin (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, fgaussOD, !Quiet);
     dchisq = ochisq - chisq;
     ochisq = chisq;
@@ -70 +76 @@
 
   for (i = 0; i < Npts; i++) {
-    ovec[0].elements[i] = fgaussOD (xvec[0].elements[i], par, Npar, dy);
+    ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, dy);
   }
   ovec[0].Nelements = Npts;
@@ -90 +96 @@
 
 /* pars: x_o, sigma, I, back */
-float fgaussOD (float x, float *par, int Npar, float *dpar) {
+opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
 
-  float z, r, f;
+  opihi_flt z, r, f;
 
   z = (x - par[0])/par[1];