IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26467


Ignore:
Timestamp:
Dec 22, 2009, 9:52:30 AM (16 years ago)
Author:
eugene
Message:

adding hermitian polynomials and functions

Location:
branches/eam_branches/20091201/Ohana/src/opihi
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c

    r26463 r26467  
    1 # include "astro.h"
     1# include "data.h"
    22
    33int hermitian1d (int argc, char **argv) {
    44 
    5   int i, order, Nvec, Poly;
    6   int *inI;
    7   float *out, inF;
    8   float mean, sigma;
     5  int i, N, order, Nvec, Poly;
     6  opihi_int *inI;
     7  opihi_flt *out, *inF;
     8  opihi_flt mean, sigma, x;
     9  double nf, norm;
    910  Vector *xvec, *yvec;
    1011
    … …  
    1516  }
    1617
    17   if (argc != 6) {
    18     gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n");
    19     return (FALSE);
    20   }
     18  if (argc != 6) goto usage;
    2119
    2220  /* select input / output buffers */
    … …  
    2725
    2826  /* gaussian parameters */
    29   mean = atof (argv[3]);
     27  mean  = atof (argv[3]);
    3028  sigma = atof (argv[4]);
    3129  order = atoi (argv[5]);
     30  if (order < 0) goto usage;
     31  if (order > 10) goto usage;
    3232
    33   out = (float *) yvec[0].elements.Flt;
    34   inF = (float *) xvec[0].elements.Flt;
    35   inI = (int *)   xvec[0].elements.Int;
     33  out = yvec[0].elements.Flt;
     34  inF = xvec[0].elements.Flt;
     35  inI = xvec[0].elements.Int;
     36
     37  nf = exp(lgamma(order + 1));
     38  norm = 1.0 / sqrt(nf*sqrt(2*M_PI)*sigma);
    3639
    3740  // a little sub-optimal : split this up with macros?
    … …  
    4447    *out = hermitian_polynomial (x, order);
    4548    if (!Poly) {
    46       *out *= exp (-0.25*x*x);
     49      *out *= norm*exp (-0.25*x*x);
    4750    }
    4851  }
    4952  return (TRUE);
     53
     54usage:
     55  gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n");
     56  gprint (GP_ERR, "  Note : only orders up to 10 are supported\n");
     57  return (FALSE);
     58
    5059}
    5160
    52   // {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\!
    53 
    54 double hermitian_polynomial (double x, int order) {
    55   double value;
    56     switch (order) {
    57       case 0:
    58         value = hermitian_00(x);
    59         break;
    60       case 1:
    61         value = hermitian_01(x);
    62         break;
    63       case 2:
    64         value = hermitian_02(x);
    65         break;
    66       case 3:
    67         value = hermitian_03(x);
    68         break;
    69       case 4:
    70         value = hermitian_04(x);
    71         break;
    72       case 5:
    73         value = hermitian_05(x);
    74         break;
    75       case 6:
    76         value = hermitian_06(x);
    77         break;
    78       case 7:
    79         value = hermitian_07(x);
    80         break;
    81       case 8:
    82         value = hermitian_08(x);
    83         break;
    84       case 9:
    85         value = hermitian_09(x);
    86         break;
    87       case 10:
    88         value = hermitian_10(x);
    89         break;
    90       default:
    91         value = NAN;
    92         break;
    93     }
    94     return value;
    95 }
    96 
    97 double hermitian_00(double x) {
    98     double value;
    99     // H_0(x) = 1
    100     value = 1;
    101     return value;
    102 }
    103 double hermitian_01(double x) {
    104     double value;
    105     // H_1(x) = x
    106     value = x;
    107     return value;
    108 }
    109 double hermitian_02(double x) {
    110     double value, x2;
    111     // H_2(x) = x^2-1
    112     x2 = x*x;
    113     value = x2 - 1.0;
    114     return value;
    115 }
    116 double hermitian_03(double x) {
    117     double value, x2;
    118     // H_3(x) = x^3-3x
    119     x2 = x*x;
    120     value = x*(x2 - 3.0);
    121     return value;
    122 }
    123 double hermitian_04(double x) {
    124     double value, x2;
    125     // H_4(x) = x^4-6x^2+3
    126     x2 = x*x;
    127     value = (x2 - 6.0)*x2 + 3.0;
    128     return value;
    129 }
    130 double hermitian_05(double x) {
    131     double value, x2;
    132     // H_5(x) = x^5-10x^3+15x
    133     x2 = x*x;
    134     value = ((x2 - 10.0)*x2 + 15.0)*x;
    135     return value;
    136 }
    137 double hermitian_06(double x) {
    138     double value, x2;
    139     // H_6(x) = x^6-15x^4+45x^2-15
    140     x2 = x*x;
    141     value = ((x2 - 15.0)*x2 + 45.0*x2) - 15.0;
    142     return value;
    143 }
    144 double hermitian_07(double x) {
    145     double value, x2;
    146     // H_7(x) = x^7-21x^5+105x^3-105x
    147     x2 = x*x;
    148     value = (((x2 - 21.0)*x2+105.0)*x2 - 105.0)*x;
    149     return value;
    150 }
    151 double hermitian_08(double x) {
    152     double value, x2;
    153     // H_8(x) = x^8-28x^6+210x^4-420x^2+105
    154     x2 = x*x;
    155     value = ((((x2 - 28.0)*x2 + 210.0)*x2 - 420.0)*x2 + 105.0);
    156     return value;
    157 }
    158 double hermitian_09(double x) {
    159     double value, x2;
    160     // H_9(x) = x^9-36x^7+378x^5-1260x^3+945x
    161     x2 = x*x;
    162     value = ((((x2 - 36.0)*x2 + 378.0)*x2 - 1260.0)*x2 + 945.0);
    163     return value;
    164 }
    165 double hermitian_10(double x) {
    166     double value, x2;
    167     // H_{10}(x) = x^{10}-45x^8+630x^6-3150x^4+4725x^2-945
    168     x2 = x*x;
    169     value = (((((x2 - 45.0)*x2 + 630.0)*x2 - 3150.0)*x2 + 4725.0)*x2 - 945.0);
    170     return value;
    171 }
     61// {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\!

  • branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c

    r26463 r26467  
    33int hermitian2d (int argc, char **argv) {
    44 
    5   int i, j, Nx, Ny, N;
     5  int i, j, N, Nx, Ny, Xorder, Yorder, Poly;
    66  float *in;
    7   double Sig_x, Sig_y, Theta;
    8   double root1, root2, R, A1, A2, A3;
    9   double Sx, Sy, Sxy;
    10   double x, y, r, f, Xo, Yo;
     7  double Xo, Yo, sigma;
     8  double x, y, xf, yf, value, norm, weight;
     9  double nf, mf;
    1110  Buffer *buf;
    1211
    13   Poly = FALSE;
     12  // optionally return the Hermitian Polynomial or the Hermitian Function
     13  Poly = FALSE;
    1414  if ((N = get_argument (argc, argv, "-poly"))) {
    1515    Poly = TRUE;
    … …  
    2626  }
    2727
    28   if (argc > 5) {
    29     gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n");
    30     return (FALSE);
    31   }
     28  if (argc != 5) goto usage;
    3229
    3330  /* select input / output buffers */
    … …  
    4340  Xorder = atof (argv[3]);
    4441  Yorder = atof (argv[4]);
     42  if (Xorder < 0) goto usage;
     43  if (Yorder < 0) goto usage;
     44  if (Xorder > 10) goto usage;
     45  if (Yorder > 10) goto usage;
     46
     47  nf = exp(lgamma(Xorder + 1));
     48  mf = exp(lgamma(Yorder + 1));
     49  norm = 1.0 / sqrt(nf*mf*2*M_PI) / sigma;
    4550
    4651  in = (float *) buf[0].matrix.buffer;
    … …  
    5863      value = xf*yf;
    5964      if (!Poly) {
    60         value *= xf*yf;
     65        weight = norm*exp(-0.25*(x*x + y*y));
     66        value *= weight;
    6167      }
    6268
    63       *in += xf*yf;
     69      *in += value;
    6470    }
    6571  }
    6672
    6773  return (TRUE);
     74
     75usage:
     76  gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n");
     77  gprint (GP_ERR, "  Note : only orders up to 10 are supported\n");
     78  return (FALSE);
    6879}

  • branches/eam_branches/20091201/Ohana/src/opihi/include/data.h

    r26418 r26467  
    165165void FreeBooks (void);
    166166
     167/* hermitian functions */
     168double hermitian_polynomial (double x, int order);
     169double hermitian_00(double x);
     170double hermitian_01(double x);
     171double hermitian_02(double x);
     172double hermitian_03(double x);
     173double hermitian_04(double x);
     174double hermitian_05(double x);
     175double hermitian_06(double x);
     176double hermitian_07(double x);
     177double hermitian_08(double x);
     178double hermitian_09(double x);
     179double hermitian_10(double x);
     180
    167181# endif

  • branches/eam_branches/20091201/Ohana/src/opihi/lib.data/Makefile

    r26417 r26467  
    2828$(SDIR)/precess.$(ARCH).o               \
    2929$(SDIR)/starfuncs.$(ARCH).o             \
     30$(SDIR)/hermitian.$(ARCH).o             \
    3031$(SDIR)/gaussian.$(ARCH).o              \
    3132$(SDIR)/graphtools.$(ARCH).o            \
Note: See TracChangeset for help on using the changeset viewer.