[Kst] extragear/graphics/kst/kst/datasources/healpix

[Ted Kisner]

SVN commit 475913 by tskisner:

readField from the healpix datasource now returns one of the 4 vectors of coordinates used in the vectorfield.  But we cannot plot vectorfields yet...

 M  +97 -14    healpix.cpp  


--- trunk/extragear/graphics/kst/kst/datasources/healpix/healpix.cpp #475912:475913
@@ -151,10 +151,10 @@
         }
         _matrixList.append(mapName);
       }
-      _fieldList.append("Vector Field Head X");
-      _fieldList.append("Vector Field Head Y");
-      _fieldList.append("Vector Field Tail X");
-      _fieldList.append("Vector Field Tail Y");
+      _fieldList.append("Vector Field Head Theta");
+      _fieldList.append("Vector Field Head Phi");
+      _fieldList.append("Vector Field Tail Theta");
+      _fieldList.append("Vector Field Tail Phi");
     } else {
       healpix_keys_free(_keys);
       healpix_strarr_free(_names, HEALPIX_FITS_MAXCOL);
@@ -187,6 +187,9 @@
   return setLastUpdateResult(KstObject::NO_CHANGE);
 }
 
+//FIXME we should really change readField so that we can compute 
+// a vector field using component maps from different files.
+// this will have to wait for the future...
 
 int HealpixSource::readField(double *v, const QString& field, int s, int n) {
   Q_UNUSED(s)
@@ -207,14 +210,13 @@
     float *mapdata;
     float *comptheta;
     float *compphi;
-    double *headsX;
-    double *headsY;
-    double *tailsX;
-    double *tailsY;
     float nullval;
     int nnull;
     size_t vecNside;
     size_t temppix;
+    int fieldnum = _fieldList.findIndex(field);
+    double theta, phi;
+    long nearest[8];
 
     // check range
     if (n <= 0) {
@@ -401,7 +403,13 @@
     // autoscale magnitude if necessary
     
     if (_autoMag) {
-      
+      double vecMag = 0.0;
+      for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
+        if (sqrt((double)(comptheta[i]*comptheta[i] + compphi[i]*compphi[i])) > vecMag) {
+          vecMag = sqrt((double)(comptheta[i]*comptheta[i] + compphi[i]*compphi[i]));
+        }
+      }
+      _maxMag = vecMag;
     }
     
     // find distance to 4 nearest neighbor pixels and 
@@ -410,9 +418,84 @@
     // This is *expensive*, but we only load vector 
     // fields infrequently, so should be ok.
     
+    // only actually compute the component of the 
+    // vectorfield that is requested.
     
+    double maxang;
+    double thetapart;
+    double phipart;
+    double headtheta;
+    double tailtheta;
+    double headphi;
+    double tailphi;
+    // only compute as many elements as requested
+    for (size_t i = 0; i < (size_t)n; i++) {
+      if (!healpix_is_fnull(comptheta[i]) && !healpix_is_fnull(compphi[i])) {
+        thetapart = (double)comptheta[i];
+        phipart = (double)compphi[i];
+        if (_mapOrder == HEALPIX_RING) {
+          healpix_pix2ang_ring(vecNside, i, &theta, &phi);
+        } else {
+          healpix_pix2ang_nest(vecNside, i, &theta, &phi);
+        }
+        //find max angle for this pixel
+        healpix_neighbors(vecNside, _mapOrder, i, nearest);
+        maxang = healpix_loc_dist(vecNside, _mapOrder, i, nearest[0]);
+        maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[1]);
+        maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[2]);
+        maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[3]);
+        maxang /= 4.0;
+        thetapart *= maxang / _maxMag;
+        phipart *= maxang / _maxMag;
+        //find coordinates of head and tail
+        headtheta = theta + 0.5 * thetapart;
+        tailtheta = theta - 0.5 * thetapart;
+        headphi = phi + 0.5 * phipart;
+        tailphi = phi - 0.5 * phipart;
+        while (headtheta > HEALPIX_PI) {
+          headtheta -= HEALPIX_PI;
+        }
+        while (headtheta < 0.0) {
+          headtheta += HEALPIX_PI;
+        }
+        while (tailtheta > HEALPIX_PI) {
+          tailtheta -= HEALPIX_PI;
+        }
+        while (tailtheta < 0.0) {
+          tailtheta += HEALPIX_PI;
+        }
+        while (headphi > 2.0 * HEALPIX_PI) {
+          headphi -= 2.0 * HEALPIX_PI;
+        }
+        while (headphi < 0.0) {
+          headphi += 2.0 * HEALPIX_PI;
+        }
+        while (tailphi > 2.0 * HEALPIX_PI) {
+          tailphi -= 2.0 * HEALPIX_PI;
+        }
+        while (tailphi < 0.0) {
+          tailphi += 2.0 * HEALPIX_PI;
+        }
+        switch (fieldnum) {
+          case 0: //head theta
+            v[i] = headtheta;
+            break;
+          case 1: //head phi
+            v[i] = headphi;
+            break;
+          case 2: //tail theta
+            v[i] = tailtheta;
+            break;
+          case 3: //tail phi
+            v[i] = tailphi;
+            break;
+          default:
+            break;
+        }
+      }
+    }
         
-    return -1;
+    return n;
   } else {
     return -1;
   }
@@ -1208,10 +1291,10 @@
   ret = healpix_fits_map_test(thealpixfile, &tNside, &tOrder, &tCoord, &tType, &tMaps);
   
   if (ret) {
-    fields.append("Vector Field Head X");
-    fields.append("Vector Field Head Y");
-    fields.append("Vector Field Tail X");
-    fields.append("Vector Field Tail Y");
+    fields.append("Vector Field Head Theta");
+    fields.append("Vector Field Head Phi");
+    fields.append("Vector Field Tail Theta");
+    fields.append("Vector Field Tail Phi");
   } else {
     return QStringList();
   }