IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 9, 2015, 3:42:26 PM (11 years ago)
Author:
eugene
Message:

move multple copies of gaussian deviates to libohana

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/src/mkcmf.c

    r39121 r39242  
    1313# define FLAGS 0x1101
    1414
    15 void gauss_init (int Nbin);
    16 double rnd_gauss (double mean, double sigma);
    1715void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1816void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    … …  
    282280  }
    283281   
    284   gauss_init (2048);
     282  gaussdev_init ();
    285283
    286284  // load test stars from a file:
    … …  
    564562
    565563
    566 static int Ngaussint = 0;
    567 static double *gaussint;
    568 
    569 extern double drand48();
    570 
    571 double gaussian (double x, double mean, double sigma) {
    572 
    573   double f;
    574 
    575   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    576 
    577   return (f);
    578 
    579 }
    580 
    581 /* integrate a gaussian from -5 sigma to +5 sigma */
    582 void gauss_init (int Nbin) {
    583  
    584   int i;
    585   double val, x, dx, dx1, dx2, dx3, df;
    586   double mean, sigma;
    587  
    588   /* no need to generate this if it already exists */
    589   if (Ngaussint == Nbin) return;
    590 
    591   // A = time(NULL);
    592   // // XXX this is expensive if called a lot (1 sec min)
    593   // // for (B = 0; A == time(NULL); B++);
    594   // B = A + 10000;
    595   // srand48(B);
    596  
    597   Ngaussint = Nbin;
    598   ALLOCATE (gaussint, double, Ngaussint + 1);
    599 
    600   val = 0;
    601   dx = 1.0 / Ngaussint;
    602   dx1 = dx / 3.0;
    603   dx2 = 2.0*dx/3.0;
    604   dx3 = dx;
    605   mean = 0.0;
    606   sigma = 1.0;
    607  
    608   for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
    609     df = (3.0*gaussian(x    , mean, sigma) +
    610           9.0*gaussian(x+dx1, mean, sigma) +
    611           9.0*gaussian(x+dx2, mean, sigma) +
    612           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    613     val += df;
    614     if (val > (i + 0.5) / (double) Ngaussint) {
    615       gaussint[i] = x + dx / 2.0;
    616       i++;
    617     }
    618   }
    619 }
    620 
    621 double rnd_gauss (double mean, double sigma) {
    622  
    623   int i;
    624   double y;
    625  
    626   y = drand48();
    627   i = Ngaussint*y;
    628   y = gaussint[i]*sigma + mean;
    629  
    630   return (y);
    631  
    632 }
    633  
    634 double int_gauss (int i) {
    635   double y;
    636   y = gaussint[i];
    637   return (y);
    638 }
    639  
    640564void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {
    641565
    … …  
    654578
    655579    if (ADDNOISE) {
    656       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    657       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    658       M[i] += fSN*rnd_gauss(0.0, 1.0);
     580      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     581      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     582      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    659583      flux = pow (10.0, -0.4*M[i]);
    660584      fSN = 1.0 / sqrt(flux);
    … …  
    700624
    701625    if (ADDNOISE) {
    702       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    703       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    704       M[i] += fSN*rnd_gauss(0.0, 1.0);
     626      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     627      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     628      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    705629      flux = pow (10.0, -0.4*M[i]);
    706630      fSN = 1.0 / sqrt(flux);
    … …  
    749673
    750674    if (ADDNOISE) {
    751       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    752       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    753       M[i] += fSN*rnd_gauss(0.0, 1.0);
     675      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     676      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     677      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    754678      flux = pow (10.0, -0.4*M[i]);
    755679      fSN = 1.0 / sqrt(flux);
    … …  
    800724
    801725    if (ADDNOISE) {
    802       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    803       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    804       M[i] += fSN*rnd_gauss(0.0, 1.0);
     726      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     727      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     728      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    805729      flux = pow (10.0, -0.4*M[i]);
    806730      fSN = 1.0 / sqrt(flux);
    … …  
    857781
    858782    if (ADDNOISE) {
    859       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    860       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    861       M[i] += fSN*rnd_gauss(0.0, 1.0);
     783      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     784      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     785      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    862786      flux = pow (10.0, -0.4*M[i]);
    863787      fSN = 1.0 / sqrt(flux);
    … …  
    918842
    919843    if (ADDNOISE) {
    920       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    921       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    922       M[i] += fSN*rnd_gauss(0.0, 1.0);
     844      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     845      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     846      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    923847      flux = pow (10.0, -0.4*M[i]);
    924848      fSN = 1.0 / sqrt(flux);
    … …  
    987911
    988912    if (ADDNOISE) {
    989       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    990       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    991       M[i] += fSN*rnd_gauss(0.0, 1.0);
     913      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     914      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     915      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    992916      flux = pow (10.0, -0.4*M[i]);
    993917      fSN = 1.0 / sqrt(flux);
    … …  
    10941018
    10951019    if (ADDNOISE) {
    1096       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1097       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1098       M[i] += fSN*rnd_gauss(0.0, 1.0);
     1020      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     1021      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     1022      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    10991023      flux = pow (10.0, -0.4*M[i]);
    11001024      fSN = 1.0 / sqrt(flux);
    … …  
    11981122
    11991123    if (ADDNOISE) {
    1200       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1201       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1202       M[i] += fSN * rnd_gauss(0.0, 1.0);
     1124      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     1125      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     1126      M[i] += fSN * gaussdev_rnd(0.0, 1.0);
    12031127      flux = pow(10.0, -0.4 * M[i]);
    12041128      fSN = 1.0 / sqrt(flux);
    … …  
    13241248
    13251249    if (ADDNOISE) {
    1326       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1327       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1328       M[i] += fSN*rnd_gauss(0.0, 1.0);
     1250      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
     1251      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
     1252      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
    13291253      flux = pow (10.0, -0.4*M[i]);
    13301254      fSN = 1.0 / sqrt(flux);
Note: See TracChangeset for help on using the changeset viewer.