IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37546


Ignore:
Timestamp:
Nov 2, 2014, 3:36:37 PM (12 years ago)
Author:
eugene
Message:

adding spherical harmonics code

Location:
branches/eam_branches/ipp-20140904/Ohana/src
Files:
5 added
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h

    r37512 r37546  
    114114  ID_MEAS_STACK_PRIMARY  = 0x00010000,  // this stack measurement is in the primary skycell
    115115  ID_MEAS_STACK_PHOT_SRC = 0x00020000,  // this measurement supplied the stack photometry
     116  ID_MEAS_ICRF_QSO       = 0x00040000,  // this measurement is an ICRF reference position
    116117} DVOMeasureFlags;
    117118
    … …  
    130131*/
    131132
     133/** these names were previously used for Average flags in old (LONEOS-era) dvo versions
     134  ID_BAD_OBJECT        = 0x00004000, // if all measurements are bad, set this bit
     135  ID_MOVING            = 0x00008000, // is a moving object
     136  ID_ROCK              = 0x0000a000, // 0x8000 + 0x2000
     137  ID_GHOST             = 0x0000c001, // 0x8000 + 0x4000 + 0x0001
     138  ID_TRAIL             = 0x0000c002, // 0x8000 + 0x4000 + 0x0002
     139  ID_BLEED             = 0x0000c003, // 0x8000 + 0x4000 + 0x0003
     140  ID_COSMIC            = 0x0000c004, // 0x8000 + 0x4000 + 0x0004
     141**/
     142
    132143/* Average.flags values -- these values are 32 bit (as of PS1_V1) */
    133144typedef enum {
    134   ID_STAR_FEW       = 0x00000001, // used within relphot: skip star
    135   ID_STAR_POOR      = 0x00000002, // used within relphot: skip star
    136   ID_PROPER         = 0x00000400, // star with large proper motion
    137   ID_TRANSIENT      = 0x00001000, // is this mutually exclusive with USNO?
    138   ID_VARIABLE       = 0x00002000, // not currently set?
    139   ID_ASTEROID       = 0x00002000, // identified with an asteroid
    140   ID_BAD_OBJECT     = 0x00004000, // if all measurements are bad, set this bit
    141   ID_MOVING         = 0x00008000, // is a moving object
    142   ID_ROCK           = 0x0000a000, // 0x8000 + 0x2000
    143   ID_GHOST          = 0x0000c001, // 0x8000 + 0x4000 + 0x0001
    144   ID_TRAIL          = 0x0000c002, // 0x8000 + 0x4000 + 0x0002
    145   ID_BLEED          = 0x0000c003, // 0x8000 + 0x4000 + 0x0003
    146   ID_COSMIC         = 0x0000c004, // 0x8000 + 0x4000 + 0x0004
    147   ID_STAR_FIT_AVE   = 0x00010000, // average position fitted
    148   ID_STAR_FIT_PM    = 0x00020000, // proper motion fitted
    149   ID_STAR_FIT_PAR   = 0x00040000, // parallax fitted
    150   ID_STAR_USE_AVE   = 0x00080000, // average position used (not PM or PAR)
    151   ID_STAR_USE_PM    = 0x00100000, // proper motion used (not AVE or PAR)
    152   ID_STAR_USE_PAR   = 0x00200000, // parallax used (not AVE or PM)
    153   ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured
    154   ID_STAR_BAD_PM    = 0x00800000, // mean astrometry could not be measured
    155   ID_OBJ_EXT        = 0x01000000, // extended in our data (eg, PS)
    156   ID_OBJ_EXT_ALT    = 0x02000000, // extended in external data (eg, 2MASS)
    157   ID_OBJ_GOOD       = 0x04000000, // good-quality measurement in our data (eg,PS)
    158   ID_OBJ_GOOD_ALT   = 0x08000000, // good-quality measurement in  external data (eg, 2MASS)
    159   ID_OBJ_GOOD_STACK = 0x10000000, // good-quality object in the stack (> 1 good stack)
     145  ID_STAR_FEW          = 0x00000001, // used within relphot: skip star
     146  ID_STAR_POOR         = 0x00000002, // used within relphot: skip star
     147  ID_ICRF_QSO          = 0x00000004, // object IDed with known ICRF quasar (may have ICRF position measurement)
     148  ID_OTHEF_QSO         = 0x00000008, // object IDed with known or suspected quasar, without ICRF reference data
     149  ID_PROPER            = 0x00000010, // star with large proper motion
     150  ID_TRANSIENT         = 0x00000020, // identified as a non-periodic (stationary) transient
     151  ID_VARIABLE          = 0x00000040, // identified as a period variable
     152  ID_ASTEROID          = 0x00000080, // identified with a known solar-system object (asteroid or other)
     153  // bits 0x00000100 - 0x00008000 are currently unused
     154  ID_STAR_FIT_AVE      = 0x00010000, // average position fitted
     155  ID_STAR_FIT_PM       = 0x00020000, // proper motion fitted
     156  ID_STAR_FIT_PAR      = 0x00040000, // parallax fitted
     157  ID_STAR_USE_AVE      = 0x00080000, // average position used (not PM or PAR)
     158  ID_STAR_USE_PM       = 0x00100000, // proper motion used (not AVE or PAR)
     159  ID_STAR_USE_PAR      = 0x00200000, // parallax used (not AVE or PM)
     160  ID_STAR_NO_ASTROM    = 0x00400000, // mean astrometry could not be measured
     161  ID_STAR_BAD_PM       = 0x00800000, // mean astrometry could not be measured
     162  ID_OBJ_EXT           = 0x01000000, // extended in our data (eg, PS)
     163  ID_OBJ_EXT_ALT       = 0x02000000, // extended in external data (eg, 2MASS)
     164  ID_OBJ_GOOD          = 0x04000000, // good-quality measurement in our data (eg,PS)
     165  ID_OBJ_GOOD_ALT      = 0x08000000, // good-quality measurement in  external data (eg, 2MASS)
     166  ID_OBJ_GOOD_STACK    = 0x10000000, // good-quality object in the stack (> 1 good stack)
    160167  ID_OBJ_SUSPECT_STACK = 0x20000000, // suspect object in the stack (> 1 good or suspect stack, < 2 good)
    161   ID_OBJ_BAD_STACK  = 0x40000000, // good-quality object in the stack (> 1 good stack)
     168  ID_OBJ_BAD_STACK     = 0x40000000, // good-quality object in the stack (> 1 good stack)
    162169} DVOAverageFlags;
    163170

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/coordops.c

    r37531 r37546  
    670670        equinox = strtod (equinoxString, &endptr);
    671671        if (endptr == equinoxString) haveEquinox = FALSE;
    672         if (!haveEquinox) {
    673           equinox = 2000.0;
    674         }
    675672      }
     673    }
     674    if (!haveEquinox) {
     675      equinox = 2000.0;
    676676    }
    677677    if (fabs (equinox - 2000.0) > 0.1) {

  • branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h

    r37538 r37546  
    448448
    449449typedef struct {
     450  double *Fr;
     451  double *Fi;
     452  int    *l;
     453  int    *m;
     454  int lmax;
     455  int Nterms;
     456} SHterms;
     457
     458SHterms *SHtermsInit (int lmax);
     459void SHtermsFree (SHterms *terms);
     460void SHtermsForRD (SHterms *terms, double R, double D);
     461void SHtermsForLM (SHterms *terms, double R, double D, int l, int m);
     462
     463typedef struct {
    450464  double *dR_B;
    451465  double *dR_E;
    452466  double *dD_B;
    453467  double *dD_E;
     468  int    *l;
     469  int    *m;
    454470  int lmax;
    455471  int Nterms;
    456472} VSHterms;
    457473
    458 VSHterms *InitVSHterms (int lmax);
    459 void getVSHterms (VSHterms *terms, double R, double D);
    460 
    461 # endif
     474VSHterms *VSHtermsInit (int lmax);
     475void VSHtermsFree (VSHterms *terms);
     476void VSHtermsForRD (VSHterms *terms, double R, double D);
     477void VSHtermsForLM (VSHterms *terms, double R, double D, int l, int m);
     478
     479# endif

  • branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/spherical.c

    r37538 r37546  
    22# include <gsl/gsl_sf_legendre.h>
    33
    4 VSHterms *InitVSHterms (int lmax) {
     4SHterms *SHtermsInit (int lmax) {
    55
    66  // 0,0  -- N = 1, S = 1
    … …  
    99  // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 (l+1)^2
    1010
     11  SHterms *terms = NULL;
     12  ALLOCATE (terms, SHterms, 1);
     13  terms->lmax = lmax;
     14  terms->Nterms = SQ(lmax+1);
     15
     16  ALLOCATE_ZERO (terms->Fr, double, terms->Nterms);
     17  ALLOCATE_ZERO (terms->Fi, double, terms->Nterms);
     18  ALLOCATE_ZERO (terms->l,     int, terms->Nterms);
     19  ALLOCATE_ZERO (terms->m,     int, terms->Nterms);
     20 
     21  int l, m;
     22  int n = 0;
     23  terms->l[n] = 0;
     24  terms->m[n] = 0;
     25  n++;
     26  for (l = 1; l <= terms->lmax; l++) {
     27    for (m = -l; m <= l; m++) {
     28      terms->l[n] = l;
     29      terms->m[n] = m;
     30      n++;
     31    }
     32  }
     33  return terms;
     34}
     35
     36void SHtermsFree (SHterms *terms) {
     37
     38  free (terms->Fr);
     39  free (terms->Fi);
     40  free (terms->l);
     41  free (terms->m);
     42
     43  free (terms);
     44  return;
     45}
     46
     47void SHtermsForRD (SHterms *terms, double R, double D) {
     48
     49  // given a point on the sky R,D
     50
     51  // for the real terms, Y_lm =
     52  // (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)]
     53  // (m <  0) : Y_{l,m} = K(l,m) P_{l,m}*sin(|m|*R)  = gsl_sphPlm(l,|m|,cos(90-D))*sin(|m|*R)
     54  // (m == 0) : Y_{l,m} = K(l,m) P_{l,m}             = gsl_sphPlm(l,|m|,cos(90-D))
     55
     56  double Rrad = R * RAD_DEG;
     57  double Drad = D * RAD_DEG;
     58
     59  int l, m;
     60
     61  int n = 0;
     62  for (l = 0; l <= terms->lmax; l++) {
     63
     64    for (m = -l; m <= l; m++) {
     65   
     66      int m_abs = abs(m);
     67      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))
     68
     69      terms->Fr[n] = aP_lm * cos (m_abs * Rrad) ;
     70      terms->Fi[n] = aP_lm * sin (m_abs * Rrad);
     71      n++;
     72    }
     73  }
     74}
     75/** I can confirm that these terms match those generated by V. Makarov's program **/
     76
     77void SHtermsForLM (SHterms *terms, double R, double D, int l, int m) {
     78
     79  double Rrad = R * RAD_DEG;
     80  double Drad = D * RAD_DEG;
     81
     82  int m_abs = abs(m);
     83  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))
     84
     85  terms->Fr[0] = aP_lm * cos (m_abs * Rrad) ;
     86  terms->Fi[0] = aP_lm * sin (m_abs * Rrad);
     87  terms->l[0] = l;
     88  terms->m[0] = m;
     89}
     90
     91VSHterms *VSHtermsInit (int lmax) {
     92
     93  // 0,0  -- N = 1, S = 1
     94  // 1,-1; 1,0; 1,1 -- N = 3, S = 4 (l+1)^2
     95  // 2,-2; 2,-1; 2,0; 2,1; 2,2 -- N = 5, S = 9 (l+1)^2
     96  // 3,-3; 3,-2; 3,-1; 3,0; 3,1; 3,2, 3,3 -- N = 7, S = 16 (l+1)^2
     97
    1198  VSHterms *terms = NULL;
    1299  ALLOCATE (terms, VSHterms, 1);
    … …  
    18105  ALLOCATE (terms->dD_B, double, terms->Nterms);
    19106  ALLOCATE (terms->dD_E, double, terms->Nterms);
    20  
     107  ALLOCATE (terms->l,    int,    terms->Nterms);
     108  ALLOCATE (terms->m,    int,    terms->Nterms);
     109 
     110  int l, m;
     111  int n = 0;
     112  terms->l[n] = 0;
     113  terms->m[n] = 0;
     114  n++;
     115  for (l = 1; l <= terms->lmax; l++) {
     116    for (m = -l; m <= l; m++) {
     117      terms->l[n] = l;
     118      terms->m[n] = m;
     119      n++;
     120    }
     121  }
     122
    21123  return terms;
    22124}
    23125
    24 void getVSHterms (VSHterms *terms, double R, double D) {
     126void VSHtermsFree (VSHterms *terms) {
     127
     128  if (!terms) return;
     129
     130  free (terms->dR_B);
     131  free (terms->dR_E);
     132  free (terms->dD_B);
     133  free (terms->dD_E);
     134  free (terms->l);
     135  free (terms->m);
     136 
     137  free (terms);
     138
     139  return;
     140}
     141
     142void VSHtermsForRD (VSHterms *terms, double R, double D) {
    25143
    26144  // given a point on the sky R,D
    … …  
    93211  n++;
    94212
    95   for (l = 1; l < terms->lmax; l++) {
     213  for (l = 1; l <= terms->lmax; l++) {
    96214
    97215    for (m = -l; m <= l; m++) {
    98216   
    99217      int m_abs = abs(m);
    100       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))
     218      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))
    101219
    102220      // note that for m == 0, S_R_lm -> 0 (due to m*)
    103221      double S_R_lm = (m > 0) ? -m * aP_lm * sin (m_abs * Rrad) : -m * aP_lm * cos (m_abs * Rrad) ;
    104222     
    105       double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad));
     223      // double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad));
     224      double aP_lp1m = gsl_sf_legendre_sphPlm(l + 1, m_abs, cos(M_PI/2 - Drad));
    106225      // reuse for aP_lm above?
    107226     
    108       double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(D);
     227      double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(Drad);
    109228      // XXX : this diverges @ D = 90.0??
     229      // this is derived from the recursion relation:
     230      // (1 - x^2) (d/dx) P_l^m(x) = (l + 1) x P_l^m(x) - (l - m + 1) P^m_{l+1}
     231      // x = cos \theta = sin D
     232      // (1 - x^2) = sin^2 \theta = cos^2 D
     233      // (d/d\theta) P_l^m(x) = [(d/dx) P_l^m(x)] [dx/d\theta]
     234      // dx/d\theta = d cos\theta / d\theta = -sin \theta = -cos D [d/dD changes the sign]
     235
     236      // (d/d\theta) P^m_l (x) = ((l + 1) sin D P_l^m(x) - (l - m + 1) P^m_{l+1}) / cos D
    110237
    111238      double S_D_lm = (m < 0) ? dPlm_dD * sin(m_abs * Rrad) : dPlm_dD * cos(m_abs * Rrad);
    … …  
    115242      terms->dR_E[n] = +S_R_lm / cos(Drad);
    116243      terms->dD_E[n] = +S_D_lm;
    117       n++;
    118     }
    119   }
    120 }
    121    
     244
     245      // 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]);
     246      // 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]);
     247      n++;
     248    }
     249  }
     250}
     251/** I can confirm that these terms match those generated by V. Makarov's program **/
     252
     253void VSHtermsForLM (VSHterms *terms, double R, double D, int l, int m) {
     254
     255  double Rrad = R * RAD_DEG;
     256  double Drad = D * RAD_DEG;
     257
     258  int m_abs = abs(m);
     259  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))
     260
     261  // note that for m == 0, S_R_lm -> 0 (due to m*)
     262  double S_R_lm = (m > 0) ? -m * aP_lm * sin (m_abs * Rrad) : -m * aP_lm * cos (m_abs * Rrad) ;
     263     
     264  double aP_lp1m = gsl_sf_legendre_Plm(l + 1, m_abs, cos(M_PI/2 - Drad));
     265  // reuse for aP_lm above?
     266     
     267  double dPlm_dD = ((l+1) * sin(Drad) * aP_lm - (l+1-m_abs)*aP_lp1m) / cos(Drad);
     268  // XXX : this diverges @ D = 90.0??
     269
     270  double S_D_lm = (m < 0) ? dPlm_dD * sin(m_abs * Rrad) : dPlm_dD * cos(m_abs * Rrad);
     271 
     272  terms->dR_B[0] = +S_D_lm;
     273  terms->dD_B[0] = -S_R_lm / cos(Drad);
     274  terms->dR_E[0] = +S_R_lm / cos(Drad);
     275  terms->dD_E[0] = +S_D_lm;
     276}

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/Makefile

    r36489 r37546  
    7171$(SRC)/star.$(ARCH).o              \
    7272$(SRC)/transform.$(ARCH).o         \
     73$(SRC)/vshimage.$(ARCH).o         \
     74$(SRC)/shimage.$(ARCH).o         \
    7375$(SRC)/imsub.$(ARCH).o             \
    7476$(SRC)/imfit.$(ARCH).o             \

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.astro/init.c

    r36489 r37546  
    6161int times                   PROTO((int, char **));
    6262int transform               PROTO((int, char **));
     63int vshimage                PROTO((int, char **));
     64int shimage                 PROTO((int, char **));
    6365
    6466static Command cmds[] = { 
    … …  
    121123  {1, "star",        star,         "star stats at rough coords"},
    122124  {1, "transform",   transform,    "geometric transformation of image"},
     125  {1, "vshimage",    vshimage,     "generate images for vector spherical harmonic terms"},
     126  {1, "shimage",     shimage,      "generate images for spherical harmonic terms"},
    123127};
    124128

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile

    r37540 r37546  
    161161$(SRC)/xsection.$(ARCH).o          \
    162162$(SRC)/vsh.$(ARCH).o       \
     163$(SRC)/vshfit.$(ARCH).o    \
     164$(SRC)/shterms.$(ARCH).o           \
     165$(SRC)/shfit.$(ARCH).o     \
    163166$(SRC)/wd.$(ARCH).o                \
    164167$(SRC)/write_vectors.$(ARCH).o     \

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c

    r37541 r37546  
    149149int vsmooth          PROTO((int, char **));
    150150int vsh              PROTO((int, char **));
     151int vshfit           PROTO((int, char **));
     152int shterms          PROTO((int, char **));
     153int shfit            PROTO((int, char **));
    151154int wd               PROTO((int, char **));
    152155int write_vectors    PROTO((int, char **));
    … …  
    317320  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    318321  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},
     322  {1, "vshfit",       vshfit,           "Vector Spherical Harmonics fits"},
     323  {1, "shterms",      shterms,          "Spherical Harmonics terms"},
     324  {1, "shfit",        shfit,            "Spherical Harmonics fits"},
    319325  {1, "wd",           wd,               "write an image to a file"},
    320326  {1, "write",        write_vectors,    "write vectors to datafile"},

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/vsh.c

    r37539 r37546  
    1919  int lmax = atoi(argv[5]);
    2020
    21   VSHterms *terms = InitVSHterms (lmax);
     21  VSHterms *terms = VSHtermsInit (lmax);
    2222
    2323  double R = atof(argv[6]);
    … …  
    2929  ResetVector (Devec, OPIHI_FLT, terms->Nterms);
    3030
    31   getVSHterms (terms, R, D);
     31  VSHtermsForRD (terms, R, D);
    3232
    3333  for (i = 0; i < terms->Nterms; i++) {
    … …  
    3838  }
    3939
    40   // FreeVSHterms (terms);
     40  VSHtermsFree (terms);
    4141
    4242  return (TRUE);
Note: See TracChangeset for help on using the changeset viewer.