Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h (revision 37546) @@ -114,4 +114,5 @@ ID_MEAS_STACK_PRIMARY = 0x00010000, // this stack measurement is in the primary skycell ID_MEAS_STACK_PHOT_SRC = 0x00020000, // this measurement supplied the stack photometry + ID_MEAS_ICRF_QSO = 0x00040000, // this measurement is an ICRF reference position } DVOMeasureFlags; @@ -130,34 +131,40 @@ */ +/** these names were previously used for Average flags in old (LONEOS-era) dvo versions + ID_BAD_OBJECT = 0x00004000, // if all measurements are bad, set this bit + ID_MOVING = 0x00008000, // is a moving object + ID_ROCK = 0x0000a000, // 0x8000 + 0x2000 + ID_GHOST = 0x0000c001, // 0x8000 + 0x4000 + 0x0001 + ID_TRAIL = 0x0000c002, // 0x8000 + 0x4000 + 0x0002 + ID_BLEED = 0x0000c003, // 0x8000 + 0x4000 + 0x0003 + ID_COSMIC = 0x0000c004, // 0x8000 + 0x4000 + 0x0004 +**/ + /* Average.flags values -- these values are 32 bit (as of PS1_V1) */ typedef enum { - ID_STAR_FEW = 0x00000001, // used within relphot: skip star - ID_STAR_POOR = 0x00000002, // used within relphot: skip star - ID_PROPER = 0x00000400, // star with large proper motion - ID_TRANSIENT = 0x00001000, // is this mutually exclusive with USNO? - ID_VARIABLE = 0x00002000, // not currently set? - ID_ASTEROID = 0x00002000, // identified with an asteroid - ID_BAD_OBJECT = 0x00004000, // if all measurements are bad, set this bit - ID_MOVING = 0x00008000, // is a moving object - ID_ROCK = 0x0000a000, // 0x8000 + 0x2000 - ID_GHOST = 0x0000c001, // 0x8000 + 0x4000 + 0x0001 - ID_TRAIL = 0x0000c002, // 0x8000 + 0x4000 + 0x0002 - ID_BLEED = 0x0000c003, // 0x8000 + 0x4000 + 0x0003 - ID_COSMIC = 0x0000c004, // 0x8000 + 0x4000 + 0x0004 - ID_STAR_FIT_AVE = 0x00010000, // average position fitted - ID_STAR_FIT_PM = 0x00020000, // proper motion fitted - ID_STAR_FIT_PAR = 0x00040000, // parallax fitted - ID_STAR_USE_AVE = 0x00080000, // average position used (not PM or PAR) - ID_STAR_USE_PM = 0x00100000, // proper motion used (not AVE or PAR) - ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) - ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured - ID_STAR_BAD_PM = 0x00800000, // mean astrometry could not be measured - ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) - ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) - ID_OBJ_GOOD = 0x04000000, // good-quality measurement in our data (eg,PS) - ID_OBJ_GOOD_ALT = 0x08000000, // good-quality measurement in external data (eg, 2MASS) - ID_OBJ_GOOD_STACK = 0x10000000, // good-quality object in the stack (> 1 good stack) + ID_STAR_FEW = 0x00000001, // used within relphot: skip star + ID_STAR_POOR = 0x00000002, // used within relphot: skip star + ID_ICRF_QSO = 0x00000004, // object IDed with known ICRF quasar (may have ICRF position measurement) + ID_OTHEF_QSO = 0x00000008, // object IDed with known or suspected quasar, without ICRF reference data + ID_PROPER = 0x00000010, // star with large proper motion + ID_TRANSIENT = 0x00000020, // identified as a non-periodic (stationary) transient + ID_VARIABLE = 0x00000040, // identified as a period variable + ID_ASTEROID = 0x00000080, // identified with a known solar-system object (asteroid or other) + // bits 0x00000100 - 0x00008000 are currently unused + ID_STAR_FIT_AVE = 0x00010000, // average position fitted + ID_STAR_FIT_PM = 0x00020000, // proper motion fitted + ID_STAR_FIT_PAR = 0x00040000, // parallax fitted + ID_STAR_USE_AVE = 0x00080000, // average position used (not PM or PAR) + ID_STAR_USE_PM = 0x00100000, // proper motion used (not AVE or PAR) + ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) + ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured + ID_STAR_BAD_PM = 0x00800000, // mean astrometry could not be measured + ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) + ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) + ID_OBJ_GOOD = 0x04000000, // good-quality measurement in our data (eg,PS) + ID_OBJ_GOOD_ALT = 0x08000000, // good-quality measurement in external data (eg, 2MASS) + ID_OBJ_GOOD_STACK = 0x10000000, // good-quality object in the stack (> 1 good stack) ID_OBJ_SUSPECT_STACK = 0x20000000, // suspect object in the stack (> 1 good or suspect stack, < 2 good) - ID_OBJ_BAD_STACK = 0x40000000, // good-quality object in the stack (> 1 good stack) + ID_OBJ_BAD_STACK = 0x40000000, // good-quality object in the stack (> 1 good stack) } DVOAverageFlags; Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/coordops.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/coordops.c (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/coordops.c (revision 37546) @@ -670,8 +670,8 @@ equinox = strtod (equinoxString, &endptr); if (endptr == equinoxString) haveEquinox = FALSE; - if (!haveEquinox) { - equinox = 2000.0; - } } + } + if (!haveEquinox) { + equinox = 2000.0; } if (fabs (equinox - 2000.0) > 0.1) { Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h (revision 37546) @@ -448,14 +448,32 @@ typedef struct { + double *Fr; + double *Fi; + int *l; + int *m; + int lmax; + int Nterms; +} SHterms; + +SHterms *SHtermsInit (int lmax); +void SHtermsFree (SHterms *terms); +void SHtermsForRD (SHterms *terms, double R, double D); +void SHtermsForLM (SHterms *terms, double R, double D, int l, int m); + +typedef struct { double *dR_B; double *dR_E; double *dD_B; double *dD_E; + int *l; + int *m; int lmax; int Nterms; } VSHterms; -VSHterms *InitVSHterms (int lmax); -void getVSHterms (VSHterms *terms, double R, double D); - -# endif +VSHterms *VSHtermsInit (int lmax); +void VSHtermsFree (VSHterms *terms); +void VSHtermsForRD (VSHterms *terms, double R, double D); +void VSHtermsForLM (VSHterms *terms, double R, double D, int l, int m); + +# endif Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c (revision 37546) @@ -2,5 +2,5 @@ # include -VSHterms *InitVSHterms (int lmax) { +SHterms *SHtermsInit (int lmax) { // 0,0 -- N = 1, S = 1 @@ -9,4 +9,91 @@ // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 (l+1)^2 + SHterms *terms = NULL; + ALLOCATE (terms, SHterms, 1); + terms->lmax = lmax; + terms->Nterms = SQ(lmax+1); + + ALLOCATE_ZERO (terms->Fr, double, terms->Nterms); + ALLOCATE_ZERO (terms->Fi, double, terms->Nterms); + ALLOCATE_ZERO (terms->l, int, terms->Nterms); + ALLOCATE_ZERO (terms->m, int, terms->Nterms); + + int l, m; + int n = 0; + terms->l[n] = 0; + terms->m[n] = 0; + n++; + for (l = 1; l <= terms->lmax; l++) { + for (m = -l; m <= l; m++) { + terms->l[n] = l; + terms->m[n] = m; + n++; + } + } + return terms; +} + +void SHtermsFree (SHterms *terms) { + + free (terms->Fr); + free (terms->Fi); + free (terms->l); + free (terms->m); + + free (terms); + return; +} + +void SHtermsForRD (SHterms *terms, double R, double D) { + + // given a point on the sky R,D + + // for the real terms, Y_lm = + // (m > 0) : Y_{l,m} = K(l,m) P_{l,m}*cos(m*R) = gsl_sphPlm(l,|m|,cos(90-D))*cos(|m|*R) [cos(x) == cos(-x)] + // (m < 0) : Y_{l,m} = K(l,m) P_{l,m}*sin(|m|*R) = gsl_sphPlm(l,|m|,cos(90-D))*sin(|m|*R) + // (m == 0) : Y_{l,m} = K(l,m) P_{l,m} = gsl_sphPlm(l,|m|,cos(90-D)) + + double Rrad = R * RAD_DEG; + double Drad = D * RAD_DEG; + + int l, m; + + int n = 0; + for (l = 0; l <= terms->lmax; l++) { + + for (m = -l; m <= l; m++) { + + int m_abs = abs(m); + double aP_lm = gsl_sf_legendre_sphPlm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D)) + + terms->Fr[n] = aP_lm * cos (m_abs * Rrad) ; + terms->Fi[n] = aP_lm * sin (m_abs * Rrad); + n++; + } + } +} +/** I can confirm that these terms match those generated by V. Makarov's program **/ + +void SHtermsForLM (SHterms *terms, double R, double D, int l, int m) { + + double Rrad = R * RAD_DEG; + double Drad = D * RAD_DEG; + + int m_abs = abs(m); + double aP_lm = gsl_sf_legendre_sphPlm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D)) + + terms->Fr[0] = aP_lm * cos (m_abs * Rrad) ; + terms->Fi[0] = aP_lm * sin (m_abs * Rrad); + terms->l[0] = l; + terms->m[0] = m; +} + +VSHterms *VSHtermsInit (int lmax) { + + // 0,0 -- N = 1, S = 1 + // 1,-1; 1,0; 1,1 -- N = 3, S = 4 (l+1)^2 + // 2,-2; 2,-1; 2,0; 2,1; 2,2 -- N = 5, S = 9 (l+1)^2 + // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 (l+1)^2 + VSHterms *terms = NULL; ALLOCATE (terms, VSHterms, 1); @@ -18,9 +105,40 @@ ALLOCATE (terms->dD_B, double, terms->Nterms); ALLOCATE (terms->dD_E, double, terms->Nterms); - + ALLOCATE (terms->l, int, terms->Nterms); + ALLOCATE (terms->m, int, terms->Nterms); + + int l, m; + int n = 0; + terms->l[n] = 0; + terms->m[n] = 0; + n++; + for (l = 1; l <= terms->lmax; l++) { + for (m = -l; m <= l; m++) { + terms->l[n] = l; + terms->m[n] = m; + n++; + } + } + return terms; } -void getVSHterms (VSHterms *terms, double R, double D) { +void VSHtermsFree (VSHterms *terms) { + + if (!terms) return; + + free (terms->dR_B); + free (terms->dR_E); + free (terms->dD_B); + free (terms->dD_E); + free (terms->l); + free (terms->m); + + free (terms); + + return; +} + +void VSHtermsForRD (VSHterms *terms, double R, double D) { // given a point on the sky R,D @@ -93,19 +211,28 @@ n++; - for (l = 1; l < terms->lmax; l++) { + for (l = 1; l <= terms->lmax; l++) { for (m = -l; m <= l; m++) { int m_abs = abs(m); - double aP_lm = gsl_sf_legendre_Plm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D)) + double aP_lm = gsl_sf_legendre_sphPlm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D)) // note that for m == 0, S_R_lm -> 0 (due to m*) double S_R_lm = (m > 0) ? -m * aP_lm * sin (m_abs * Rrad) : -m * aP_lm * cos (m_abs * Rrad) ; - double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad)); + // double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad)); + double aP_lp1m = gsl_sf_legendre_sphPlm(l + 1, m_abs, cos(M_PI/2 - Drad)); // reuse for aP_lm above? - double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(D); + double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(Drad); // XXX : this diverges @ D = 90.0?? + // this is derived from the recursion relation: + // (1 - x^2) (d/dx) P_l^m(x) = (l + 1) x P_l^m(x) - (l - m + 1) P^m_{l+1} + // x = cos \theta = sin D + // (1 - x^2) = sin^2 \theta = cos^2 D + // (d/d\theta) P_l^m(x) = [(d/dx) P_l^m(x)] [dx/d\theta] + // dx/d\theta = d cos\theta / d\theta = -sin \theta = -cos D [d/dD changes the sign] + + // (d/d\theta) P^m_l (x) = ((l + 1) sin D P_l^m(x) - (l - m + 1) P^m_{l+1}) / cos D double S_D_lm = (m < 0) ? dPlm_dD * sin(m_abs * Rrad) : dPlm_dD * cos(m_abs * Rrad); @@ -115,7 +242,35 @@ terms->dR_E[n] = +S_R_lm / cos(Drad); terms->dD_E[n] = +S_D_lm; - n++; - } - } -} - + + // fprintf (stderr, "%d : %d,%d : %f %f : %f %f %f %f\n", n, l, m, S_R_lm, S_D_lm, terms->dR_B[n], terms->dR_E[n], terms->dD_B[n], terms->dD_E[n]); + // fprintf (stderr, "%d : %d,%d : %f %f %f %f\n", n-1, l, m, terms->dR_B[n], terms->dR_E[n], terms->dD_B[n], terms->dD_E[n]); + n++; + } + } +} +/** I can confirm that these terms match those generated by V. Makarov's program **/ + +void VSHtermsForLM (VSHterms *terms, double R, double D, int l, int m) { + + double Rrad = R * RAD_DEG; + double Drad = D * RAD_DEG; + + int m_abs = abs(m); + double aP_lm = gsl_sf_legendre_Plm(l, m_abs, cos(M_PI/2 - Drad)); // alt-P_lm = P_lm (l, |m|, cos(90 - D)) + + // note that for m == 0, S_R_lm -> 0 (due to m*) + double S_R_lm = (m > 0) ? -m * aP_lm * sin (m_abs * Rrad) : -m * aP_lm * cos (m_abs * Rrad) ; + + double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad)); + // reuse for aP_lm above? + + double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(Drad); + // XXX : this diverges @ D = 90.0?? + + double S_D_lm = (m < 0) ? dPlm_dD * sin(m_abs * Rrad) : dPlm_dD * cos(m_abs * Rrad); + + terms->dR_B[0] = +S_D_lm; + terms->dD_B[0] = -S_R_lm / cos(Drad); + terms->dR_E[0] = +S_R_lm / cos(Drad); + terms->dD_E[0] = +S_D_lm; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/Makefile (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/Makefile (revision 37546) @@ -71,4 +71,6 @@ $(SRC)/star.$(ARCH).o \ $(SRC)/transform.$(ARCH).o \ +$(SRC)/vshimage.$(ARCH).o \ +$(SRC)/shimage.$(ARCH).o \ $(SRC)/imsub.$(ARCH).o \ $(SRC)/imfit.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/init.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/init.c (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/init.c (revision 37546) @@ -61,4 +61,6 @@ int times PROTO((int, char **)); int transform PROTO((int, char **)); +int vshimage PROTO((int, char **)); +int shimage PROTO((int, char **)); static Command cmds[] = { @@ -121,4 +123,6 @@ {1, "star", star, "star stats at rough coords"}, {1, "transform", transform, "geometric transformation of image"}, + {1, "vshimage", vshimage, "generate images for vector spherical harmonic terms"}, + {1, "shimage", shimage, "generate images for spherical harmonic terms"}, }; Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/shimage.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/shimage.c (revision 37546) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/shimage.c (revision 37546) @@ -0,0 +1,64 @@ +# include "data.h" + +int shimage (int argc, char **argv) { + + int ix, iy; + Buffer *FrIm, *FiIm, *src; + + if (argc != 6) { + gprint (GP_ERR, "USAGE: shimage (image) l m Fr Fi\n"); + gprint (GP_ERR, " set the Fr and Fi values for pixels in the given image for (l,m)\n"); + return (FALSE); + } + + if ((src = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + Coords coords; + GetCoords (&coords, &src[0].header); + if (!strcmp(&coords.ctype[4], "-WRP")) { + fprintf (stderr, "mosaic astrometry not yet supported for vshimage\n"); + return FALSE; + } + + int l = atoi(argv[2]); + int m = atoi(argv[3]); + + if ((FrIm = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((FiIm = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + int Nx = src[0].header.Naxis[0]; + int Ny = src[0].header.Naxis[1]; + + gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0); + + // l=0 allocates space for a single value + SHterms *terms = SHtermsInit (0); + + float *Fr = (float *) FrIm[0].matrix.buffer; + float *Fi = (float *) FiIm[0].matrix.buffer; + + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + + int Npix = ix + iy*Nx; + + double R, D; + int status = XY_to_RD (&R, &D, ix, iy, &coords); + + if (!status) { + Fr[Npix] = NAN; + Fi[Npix] = NAN; + continue; + } + + SHtermsForLM (terms, R, D, l, m); + + Fr[Npix] = terms->Fr[0]; + Fi[Npix] = terms->Fi[0]; + } + } + SHtermsFree (terms); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/vshimage.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/vshimage.c (revision 37546) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/vshimage.c (revision 37546) @@ -0,0 +1,74 @@ +# include "data.h" + +int vshimage (int argc, char **argv) { + + int ix, iy; + Buffer *Rbim, *Reim, *Dbim, *Deim, *src; + + if (argc != 8) { + gprint (GP_ERR, "USAGE: vsh image l m dRb dRe dDb dDe\n"); + gprint (GP_ERR, " set the dRb, dRe, dDb, dDe values for pixels in the given image for (l,m)\n"); + return (FALSE); + } + + if ((src = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + Coords coords; + GetCoords (&coords, &src[0].header); + if (!strcmp(&coords.ctype[4], "-WRP")) { + fprintf (stderr, "mosaic astrometry not yet supported for vshimage\n"); + return FALSE; + } + + int l = atoi(argv[2]); + int m = atoi(argv[3]); + + if ((Rbim = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((Reim = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((Dbim = SelectBuffer (argv[6], ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((Deim = SelectBuffer (argv[7], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + int Nx = src[0].header.Naxis[0]; + int Ny = src[0].header.Naxis[1]; + + gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0); + + // l=0 allocates space for a single value + VSHterms *terms = VSHtermsInit (0); + + float *Rb = (float *) Rbim[0].matrix.buffer; + float *Re = (float *) Reim[0].matrix.buffer; + float *Db = (float *) Dbim[0].matrix.buffer; + float *De = (float *) Deim[0].matrix.buffer; + + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + + int Npix = ix + iy*Nx; + + double R, D; + int status = XY_to_RD (&R, &D, ix, iy, &coords); + + if (!status) { + Rb[Npix] = NAN; + Re[Npix] = NAN; + Db[Npix] = NAN; + De[Npix] = NAN; + continue; + } + + VSHtermsForLM (terms, R, D, l, m); + + Rb[Npix] = terms->dR_B[0]; + Re[Npix] = terms->dR_E[0]; + Db[Npix] = terms->dD_B[0]; + De[Npix] = terms->dD_E[0]; + } + } + VSHtermsFree (terms); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile (revision 37546) @@ -161,4 +161,7 @@ $(SRC)/xsection.$(ARCH).o \ $(SRC)/vsh.$(ARCH).o \ +$(SRC)/vshfit.$(ARCH).o \ +$(SRC)/shterms.$(ARCH).o \ +$(SRC)/shfit.$(ARCH).o \ $(SRC)/wd.$(ARCH).o \ $(SRC)/write_vectors.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c (revision 37546) @@ -149,4 +149,7 @@ int vsmooth PROTO((int, char **)); int vsh PROTO((int, char **)); +int vshfit PROTO((int, char **)); +int shterms PROTO((int, char **)); +int shfit PROTO((int, char **)); int wd PROTO((int, char **)); int write_vectors PROTO((int, char **)); @@ -317,4 +320,7 @@ {1, "vzload", vzload, "load vectors as overlay on image display (scaled points)"}, {1, "vsh", vsh, "Vector Spherical Harmonics"}, + {1, "vshfit", vshfit, "Vector Spherical Harmonics fits"}, + {1, "shterms", shterms, "Spherical Harmonics terms"}, + {1, "shfit", shfit, "Spherical Harmonics fits"}, {1, "wd", wd, "write an image to a file"}, {1, "write", write_vectors, "write vectors to datafile"}, Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c (revision 37546) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c (revision 37546) @@ -0,0 +1,69 @@ +# include "data.h" + +int shfit (int argc, char **argv) { + + int i, j; + Vector *Rvec, *Dvec, *Fvec, *Lvec, *Mvec, *Frvec, *Fivec; + double *Fr, *Fi; + + if (argc != 9) { + gprint (GP_ERR, "USAGE: shfit R D value lmax Lout Mout Fr Fi\n"); + gprint (GP_ERR, " fit the scalar field value (at points R,D) to spherical harmonics up to the given lmax\n"); + return (FALSE); + } + + if ((Rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((Dvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((Fvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + int lmax = atoi(argv[4]); + + if ((Lvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((Mvec = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((Frvec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((Fivec = SelectVector (argv[8], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + SHterms *terms = SHtermsInit (lmax); + + ALLOCATE_ZERO (Fr, double, terms->Nterms); + ALLOCATE_ZERO (Fi, double, terms->Nterms); + + ResetVector (Frvec, OPIHI_FLT, terms->Nterms); + ResetVector (Fivec, OPIHI_FLT, terms->Nterms); + ResetVector (Lvec, OPIHI_FLT, terms->Nterms); + ResetVector (Mvec, OPIHI_FLT, terms->Nterms); + + // measure the dot product \sum(F_i * Ylm_i) + for (i = 0; i < Rvec->Nelements; i++) { + + SHtermsForRD (terms, Rvec->elements.Flt[i], Dvec->elements.Flt[i]); + + double Fv = Fvec->elements.Flt[i]; + + for (j = 0; j < terms->Nterms; j++) { + Fr[j] += Fv * terms->Fr[j]; + Fi[j] += Fv * terms->Fi[j]; + } + } + + // for (j = 0; j < terms->Nterms; j++) { + // fprintf (stderr, "%d : %d %d : %f %f\n", j, terms->l[j], terms->m[j], Fr[j], Fi[j]); + // } + + for (j = 0; j < terms->Nterms; j++) { + Fr[j] /= Fvec->Nelements; + Fi[j] /= Fvec->Nelements; + } + + for (j = 0; j < terms->Nterms; j++) { + // fprintf (stderr, "%d : %d %d : %f %f\n", j, terms->l[j], terms->m[j], 4*M_PI*Fr[j], 4*M_PI*Fi[j]); + Lvec[0].elements.Flt[j] = terms->l[j]; + Mvec[0].elements.Flt[j] = terms->m[j]; + Frvec[0].elements.Flt[j] = 4*M_PI*Fr[j]; + Fivec[0].elements.Flt[j] = 4*M_PI*Fi[j]; + } + + SHtermsFree (terms); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shterms.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shterms.c (revision 37546) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shterms.c (revision 37546) @@ -0,0 +1,43 @@ +# include "data.h" + +int shterms (int argc, char **argv) { + + int i; + Vector *Frvec, *Fivec, *lvec, *mvec; + + if (argc != 8) { + gprint (GP_ERR, "USAGE: shterms Fr Fi l m lmax R D\n"); + gprint (GP_ERR, " set the Fr and Fi values for sh terms up to lmax at R,D\n"); + return (FALSE); + } + + if ((Frvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((Fivec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((lvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((mvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + int lmax = atoi(argv[5]); + + SHterms *terms = SHtermsInit (lmax); + + double R = atof(argv[6]); + double D = atof(argv[7]); + + ResetVector (Frvec, OPIHI_FLT, terms->Nterms); + ResetVector (Fivec, OPIHI_FLT, terms->Nterms); + ResetVector (lvec, OPIHI_FLT, terms->Nterms); + ResetVector (mvec, OPIHI_FLT, terms->Nterms); + + SHtermsForRD (terms, R, D); + + for (i = 0; i < terms->Nterms; i++) { + Frvec[0].elements.Flt[i] = terms->Fr[i]; + Fivec[0].elements.Flt[i] = terms->Fi[i]; + lvec[0].elements.Flt[i] = terms->l[i]; + mvec[0].elements.Flt[i] = terms->m[i]; + } + + SHtermsFree (terms); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vsh.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vsh.c (revision 37545) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vsh.c (revision 37546) @@ -19,5 +19,5 @@ int lmax = atoi(argv[5]); - VSHterms *terms = InitVSHterms (lmax); + VSHterms *terms = VSHtermsInit (lmax); double R = atof(argv[6]); @@ -29,5 +29,5 @@ ResetVector (Devec, OPIHI_FLT, terms->Nterms); - getVSHterms (terms, R, D); + VSHtermsForRD (terms, R, D); for (i = 0; i < terms->Nterms; i++) { @@ -38,5 +38,5 @@ } - // FreeVSHterms (terms); + VSHtermsFree (terms); return (TRUE); Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vshfit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vshfit.c (revision 37546) +++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vshfit.c (revision 37546) @@ -0,0 +1,62 @@ +# include "data.h" + +int vshfit (int argc, char **argv) { + + int i, j; + Vector *Rvec, *Dvec, *dRvec, *dDvec; + double *Re, *Rb, *De, *Db; + + if (argc != 6) { + gprint (GP_ERR, "USAGE: vshfit Rvec Dvec dRvec dDvec lmax\n"); + gprint (GP_ERR, " fit the vector field dRvec, dDvec (at points Rvec,Dvec) to vector spherical harmonics up to the given lmax\n"); + return (FALSE); + } + + if ((Rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((Dvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((dRvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((dDvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + int lmax = atoi(argv[5]); + + VSHterms *terms = VSHtermsInit (lmax); + + ALLOCATE_ZERO (Re, double, terms->Nterms); + ALLOCATE_ZERO (Rb, double, terms->Nterms); + ALLOCATE_ZERO (De, double, terms->Nterms); + ALLOCATE_ZERO (Db, double, terms->Nterms); + + for (i = 0; i < Rvec->Nelements; i++) { + + VSHtermsForRD (terms, Rvec->elements.Flt[i], Dvec->elements.Flt[i]); + + double dR = dRvec->elements.Flt[i]; + double dD = dDvec->elements.Flt[i]; + + for (j = 0; j < terms->Nterms; j++) { + Rb[j] += dR * terms->dR_B[j]; + Re[j] += dR * terms->dR_E[j]; + Db[j] += dD * terms->dD_B[j]; + De[j] += dD * terms->dD_E[j]; + } + } + + for (j = 0; j < terms->Nterms; j++) { + fprintf (stderr, "%d : %d %d : %f %f %f %f\n", j, terms->l[j], terms->m[j], Rb[j], Re[j], Db[j], De[j]); + } + + for (j = 0; j < terms->Nterms; j++) { + Rb[j] /= terms->Nterms; + Re[j] /= terms->Nterms; + Db[j] /= terms->Nterms; + De[j] /= terms->Nterms; + } + + for (j = 0; j < terms->Nterms; j++) { + fprintf (stderr, "%d : %d %d : %f %f %f %f\n", j, terms->l[j], terms->m[j], Rb[j], Re[j], Db[j], De[j]); + } + + VSHtermsFree (terms); + + return (TRUE); +}