Changeset 32978

Timestamp:

Dec 20, 2011, 8:13:23 AM (15 years ago)

Author:

eugene

Message:

mods to spex1dgas

File:

• branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex1dgas.c (modified) (6 diffs)

branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex1dgas.c

@@ -113 +113 @@
   }
  
-  if (argc != 9) goto usage;
-
+  if (argc != 11) goto usage;
+
+  // XXX enforce matching lengths on the three vectors
   if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
   if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
   if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+
   int Niter = atoi (argv[4]);
   int nCloseMax = atoi (argv[5]);
@@ -123 +125 @@
   float farFrac = atof (argv[7]);
   float maxPressure = atof (argv[8]);
-  // XXX enforce matching lengths on the three vectors
+
+  int idx1 = atoi(argv[9]);
+  int idx2 = atoi(argv[10]);
 
   CastVector (index1, OPIHI_INT);
@@ -214 +218 @@
   
   // need to get pretty close to a good starting point.  start with a guess using the first entry & its extremum
-  int idx1 = 0;
-  int idx2 = object[idx1].index[object[idx1].Nindex-1];
-  float A = object[idx1].Dtgt[object[idx1].Nindex-1];
+  // int idx1 = object[0].index[(int)(0.75*object[0].Nindex)];
+  // int idx2 = object[idx1].index[(int)(0.75*object[idx1].Nindex)];
+  float A = NAN;
+  for (i = 0; isnan (A) && (i < object[idx1].Nindex); i++) {
+      if (object[idx1].index[i] == idx2) {
+          A = object[idx1].Dtgt[i];
+      }
+  }
+  fprintf (stderr, "idx1: %d, idx2: %d, A: %f\n", idx1, idx2, A);
   object[idx1].Xo = 0.0;
   object[idx1].Yo = 0.0;
@@ -227 +237 @@
   object[idx2].dPdY = 0.0;
                  
+  Vector *tmpB = SelectVector ("sp1d_B", ANYVECTOR, TRUE); 
+  Vector *tmpC = SelectVector ("sp1d_C", ANYVECTOR, TRUE); 
+  ResetVector (tmpB, OPIHI_FLT, Nobject);
+  ResetVector (tmpC, OPIHI_FLT, Nobject);
+
   for (i = 0; i < Nobject; i++) {
     if (i == idx1) continue;
@@ -238 +253 @@
     if (isnan(B) || isnan(C)) abort();
 
+    // C^2 = A^2 + B^2 - 2AB cos(t)
+    // Xo = B cos(t)
+    // Xo = A^2 + B^2 - C^2 / 2 A
     float Xo = (SQ(A) + SQ(B) - SQ(C)) / (2*A);
 
     object[i].Xo = Xo;
     object[i].dPdX = 0.0;
-  }
+
+    tmpB->elements.Flt[i] = B;
+    tmpC->elements.Flt[i] = C;
+  }
+  tmpB->elements.Flt[idx1] = 0;
+  tmpB->elements.Flt[idx2] = A;
+  tmpC->elements.Flt[idx1] = A;
+  tmpC->elements.Flt[idx2] = 0;
 
   if (Niter >= -1) {
@@ -389 +414 @@
   
 usage:
-  gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure)\n");
+  gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure) (idx1) (idx2)\n");
   return FALSE;
 }