Changeset 3349 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.cc


Ignore:
Timestamp:
Oct 11, 2007, 4:39:58 PM (18 years ago)
Author:
cmv
Message:
  • gros changements dans la structure de la classe GeneFluct3D (constructeur et logique d'aloocation memoire, init_fftw etc...)
  • suppression des valeurs de masse<0 mises a -999. directement mises a zero
  • suppression de TurnMass2HIMass qui fait maintenant la meme chose que TurnMass2MeanNumber
  • legere restructuration de cmvobserv3d.cc pour compat.

cmv 11/10/2007

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3331 r3349  
    2424#include "genefluct3d.h"
    2525
    26 #define FFTW_THREAD
     26#define WITH_FFTW_THREAD
    2727
    2828#define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
    … …  
    3131
    3232//-------------------------------------------------------
    33 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
    34   : Nx_(0) , Ny_(0) , Nz_(0)
    35   , lp_(0)
    36   , array_allocated_(false) , T_(T)
    37   , cosmo_(NULL) , growth_(NULL)
    38   , redsh_ref_(-999.), kredsh_ref_(0.), dred_ref_(-999.)
    39   , loscom_ref_(-999.), dtrc_ref_(-999.), dlum_ref_(-999.), dang_ref_(-999.)
    40   , nu_ref_(-999.), dnu_ref_ (-999.)
    41   , loscom_min_(-999.), loscom_max_(-999.)
    42   , loscom2zred_min_(0.) , loscom2zred_max_(0.)
    43 {
     33GeneFluct3D::GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz
     34                        ,unsigned short nthread,int lp)
     35{
     36 init_default();
     37
     38 lp_ = lp;
     39 nthread_ = nthread;
     40
     41 setsize(nx,ny,nz,dx,dy,dz);
     42 setalloc();
     43 setpointers(false);
     44 init_fftw();
     45}
     46
     47GeneFluct3D::GeneFluct3D(unsigned short nthread)
     48{
     49 init_default();
     50 setsize(2,2,2,1.,1.,1.);
     51 nthread_ = nthread;
     52 setalloc();
     53 setpointers(false);
     54 init_fftw();
     55}
     56
     57GeneFluct3D::~GeneFluct3D(void)
     58{
     59 delete_fftw();
     60}
     61
     62//-------------------------------------------------------
     63void GeneFluct3D::init_default(void)
     64{
     65 Nx_ = Ny_ = Nz_ = 0;
     66 is_set_fftw_plan = false;
     67 nthread_ = 0;
     68 lp_ = 0;
     69 array_allocated_ = false;
     70 cosmo_ = NULL;
     71 growth_ = NULL;
     72 redsh_ref_ = -999.;
     73 kredsh_ref_ = 0.;
     74 dred_ref_ = -999.;
     75 loscom_ref_ = -999.;
     76 dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
     77 nu_ref_ = dnu_ref_ = -999.;
     78 loscom_min_ = loscom_max_ = -999.;
     79 loscom2zred_min_ = loscom2zred_max_ = 0.;
    4480 xobs_[0] = xobs_[1] = xobs_[2] = 0.;
    4581 zred_.resize(0);
    4682 loscom_.resize(0);
    4783 loscom2zred_.resize(0);
    48  SetNThread();
    49 }
    50 
    51 GeneFluct3D::~GeneFluct3D(void)
    52 {
    53  fftw_destroy_plan(pf_);
    54  fftw_destroy_plan(pb_);
    55 #ifdef FFTW_THREAD
    56  if(nthread_>0) fftw_cleanup_threads();
    57 #endif
    58 }
    59 
    60 //-------------------------------------------------------
    61 void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
    62 {
    63   setsize(nx,ny,nz,dx,dy,dz);
    64   setalloc();
    65   setpointers(false);
    66   init_fftw();
    67 }
    68 
    69 void GeneFluct3D::SetObservator(double redshref,double kredshref)
    70 // L'observateur est au redshift z=0
    71 //               est situe sur la "perpendiculaire" a la face x,y
    72 //                         issue du centre de cette face
    73 // Il faut positionner le cube sur l'axe des z cad des redshifts:
    74 //     redshref  = redshift de reference
    75 //                 Si redshref<0 alors redshref=0
    76 //     kredshref = indice (en double) correspondant a ce redshift
    77 //                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
    78 // Exemple: redshref=1.5 kredshref=250.75
    79 //    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
    80 {
    81  if(redshref<0.) redshref = 0.;
    82  if(kredshref<0.) {
    83    if(Nz_<=0) {
    84      char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 SetSize should be called first";
    85      cout<<bla<<endl; throw ParmError(bla);
    86    }
    87    kredshref = Nz_/2.;
    88  }
    89  redsh_ref_  = redshref;
    90  kredsh_ref_ = kredshref;
    91  if(lp_>0)
    92    cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
    93 }
    94 
    95 void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
    96 {
    97  cosmo_ = &cosmo;
    98  if(lp_>1) cosmo_->Print();
    99 }
    100 
    101 void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
    102 {
    103  growth_ = &growth;
    10484}
    10585
    … …  
    171151}
    172152
    173 void GeneFluct3D::check_array_alloc(void)
    174 // Pour tester si le tableau T_ est alloue
    175 {
    176  if(array_allocated_) return;
    177  char bla[90];
    178  sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
    179  cout<<bla<<endl; throw ParmError(bla);
    180 }
    181 
    182153void GeneFluct3D::init_fftw(void)
    183154{
     155 if( is_set_fftw_plan ) delete_fftw();
     156
    184157 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
    185158 if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
    186 #ifdef FFTW_THREAD
     159#ifdef WITH_FFTW_THREAD
    187160 if(nthread_>0) {
    188161   cout<<"...Computing with "<<nthread_<<" threads"<<endl;
    … …  
    195168 if(lp_>1) cout<<"...backward plan"<<endl;
    196169 pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
     170 is_set_fftw_plan = true;
     171}
     172
     173void GeneFluct3D::delete_fftw(void)
     174{
     175 if( !is_set_fftw_plan ) return;
     176 fftw_destroy_plan(pf_);
     177 fftw_destroy_plan(pb_);
     178#ifdef WITH_FFTW_THREAD
     179 if(nthread_>0) fftw_cleanup_threads();
     180#endif
     181 is_set_fftw_plan = false;
     182}
     183
     184void GeneFluct3D::check_array_alloc(void)
     185// Pour tester si le tableau T_ est alloue
     186{
     187 if(array_allocated_) return;
     188 char bla[90];
     189 sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
     190 cout<<bla<<endl; throw ParmError(bla);
    197191}
    198192
    199193//-------------------------------------------------------
     194void GeneFluct3D::SetObservator(double redshref,double kredshref)
     195// L'observateur est au redshift z=0
     196//               est situe sur la "perpendiculaire" a la face x,y
     197//                         issue du centre de cette face
     198// Il faut positionner le cube sur l'axe des z cad des redshifts:
     199//     redshref  = redshift de reference
     200//                 Si redshref<0 alors redshref=0
     201//     kredshref = indice (en double) correspondant a ce redshift
     202//                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
     203// Exemple: redshref=1.5 kredshref=250.75
     204//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
     205{
     206 if(redshref<0.) redshref = 0.;
     207 if(kredshref<0.) {
     208   if(Nz_<=0) {
     209     char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 define cube geometry first";
     210     cout<<bla<<endl; throw ParmError(bla);
     211   }
     212   kredshref = Nz_/2.;
     213 }
     214 redsh_ref_  = redshref;
     215 kredsh_ref_ = kredshref;
     216 if(lp_>0)
     217   cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
     218}
     219
     220void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
     221{
     222 cosmo_ = &cosmo;
     223 if(lp_>1) cosmo_->Print();
     224}
     225
     226void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
     227{
     228 growth_ = &growth;
     229}
     230
    200231long GeneFluct3D::LosComRedshift(double zinc,long npoints)
    201232// Given a position of the cube relative to the observer
    … …  
    11641195}
    11651196
    1166 double  GeneFluct3D::TurnMass2HIMass(double m_by_mpc3)
    1167 {
    1168  if(lp_>0) cout<<"--- TurnMass2HIMass : "<<m_by_mpc3<<" Msol/Mpc^3"<<endl;
    1169 
    1170  double mall=0., mgood=0.;
    1171  int_8  nall=0,  ngood=0;
    1172  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    1173    int_8 ip = IndexR(i,j,l);
    1174    mall += data_[ip]; nall++;
    1175    if(data_[ip]>0.) {mgood += data_[ip]; ngood++;}
    1176  }
    1177  if(ngood>0) mgood /= (double)ngood;
    1178  if(nall>0)  mall  /= (double)nall;
    1179  if(lp_>0) cout<<"...ngood="<<ngood<<" mgood="<<mgood
    1180                <<",  nall="<<nall<<" mall="<<mall<<endl;
    1181  if(ngood<=0 || mall<=0.) {
    1182    cout<<"TurnMass2HIMass_Error: ngood="<<ngood<<" <=0 || mall="<<mall<<" <=0"<<endl;
    1183    throw RangeCheckError("TurnMass2HIMass_Error: ngood<=0 || mall<=0");
    1184  }
    1185 
    1186  // On doit mettre m*Vol masse de HI dans notre survey
    1187  // On en met uniquement dans les pixels de masse >0.
    1188  // On MET a zero les pixels <0
    1189  // On renormalise sur les pixels>0 pour qu'on ait m*Vol masse de HI
    1190  //   comme on ne prend que les pixels >0, on doit normaliser
    1191  //   a la moyenne de <1+d_rho/rho> sur ces pixels
    1192  //   (rappel sur tout les pixels <1+d_rho/rho>=1)
    1193  // masse de HI a mettre ds 1 px:
    1194  double dm = m_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
    1195  if(lp_>0) cout<<"...HI mass density move from "
    1196                <<m_by_mpc3*Vol_/double(NRtot_)<<" to "<<dm<<" / pixel"<<endl;
    1197 
    1198  double sum = 0.;
    1199  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    1200    int_8 ip = IndexR(i,j,l);
    1201    if(data_[ip]<=0.) data_[ip] = 0.;
    1202    else {
    1203      data_[ip] *= dm; // par coherence on multiplie aussi les <=0
    1204      sum += data_[ip];  // mais on ne les compte pas
    1205    }
    1206  }
    1207 
    1208  if(lp_>0) cout<<sum<<"...mass HI put into survey / "<<m_by_mpc3*Vol_<<endl;
    1209 
    1210  return sum;
    1211 }
    1212 
    1213 double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
    1214 // do NOT treate negative or nul values
    1215 {
    1216  if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<n_by_mpc3<<" gal/Mpc^3"<<endl;
     1197double GeneFluct3D::TurnMass2MeanNumber(double val_by_mpc3)
     1198{
     1199 if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;
    12171200
    12181201 double mall=0., mgood=0.;
    … …  
    12321215 }
    12331216
    1234  // On doit mettre n*Vol galaxies dans notre survey
     1217 // On doit mettre n*Vol galaxies (ou Msol) dans notre survey
    12351218 // On en met uniquement dans les pixels de masse >0.
    1236  // On NE met PAS a zero les pixels <0
    1237  // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
     1219 // On MET a zero les pixels <0
     1220 // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies (ou Msol)
    12381221 //   comme on ne prend que les pixels >0, on doit normaliser
    12391222 //   a la moyenne de <1+d_rho/rho> sur ces pixels
    12401223 //   (rappel sur tout les pixels <1+d_rho/rho>=1)
    1241  // nb de gal a mettre ds 1 px:
    1242  double dn = n_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
    1243  if(lp_>0) cout<<"...galaxy density move from "
    1244                <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
     1224 // nb de gal (ou Msol) a mettre ds 1 px:
     1225 double dn = val_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
     1226 if(lp_>0) cout<<"...quantity density move from "
     1227               <<val_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
    12451228
    12461229 double sum = 0.;
    12471230 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    12481231   int_8 ip = IndexR(i,j,l);
    1249    data_[ip] *= dn; // par coherence on multiplie aussi les <=0
    1250    if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
    1251  }
    1252 
    1253  if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
     1232   if(data_[ip]<=0.) data_[ip] = 0.;
     1233   else {
     1234     data_[ip] *= dn;
     1235     sum += data_[ip];
     1236   }
     1237 }
     1238
     1239 if(lp_>0) cout<<sum<<"...quantity put into survey / "<<val_by_mpc3*Vol_<<endl;
    12541240
    12551241 return sum;
    … …  
    13291315}
    13301316
     1317void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
     1318// add noise to every pixels (meme les <=0 !)
     1319// type_evol = 0 : pas d'evolution de la puissance du bruit
     1320//             1 : evolution de la puissance du bruit avec la distance a l'observateur
     1321//             2 : evolution de la puissance du bruit avec la distance du plan Z
     1322//                 (tous les plans Z sont mis au meme redshift z de leur milieu)
     1323{
     1324 if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
     1325 check_array_alloc();
     1326
     1327 if(type_evol<0) type_evol = 0;
     1328 if(type_evol>2) {
     1329   char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
     1330   cout<<bla<<endl; throw ParmError(bla);
     1331 }
     1332
     1333 vector<double> correction;
     1334 InterpFunc *intercor = NULL;
     1335
     1336 if(type_evol>0) {
     1337   // Sigma_Noise(en mass) :
     1338   //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
     1339   //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
     1340   //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
     1341   //                                    nlobe ~ 1/Dtrcom^2
     1342   //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
     1343   if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
     1344     char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
     1345     cout<<bla<<endl; throw ParmError(bla);
     1346   }
     1347   InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
     1348   long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
     1349   for(long i=0;i<nsz;i++) {
     1350     double d = interpinv.X(i);
     1351     double zred = interpinv(d);
     1352     double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
     1353     double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
     1354     double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
     1355     double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
     1356     double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
     1357     if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
     1358       cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
     1359           <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
     1360     correction.push_back(corr);
     1361   }
     1362   intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
     1363 }
     1364
     1365 double corrlim[2] = {1.,1.};
     1366 for(long i=0;i<Nx_;i++) {
     1367   double dx2 = DXcom(i); dx2 *= dx2;
     1368   for(long j=0;j<Ny_;j++) {
     1369     double dy2 = DYcom(j); dy2 *= dy2;
     1370     for(long l=0;l<Nz_;l++) {
     1371       double corr = 1.;
     1372       if(type_evol>0) {
     1373         double dz = DZcom(l);
     1374         if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
     1375           else dz = fabs(dz); // tous les plans Z au meme redshift
     1376         corr = (*intercor)(dz);
     1377         if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
     1378       }
     1379       int_8 ip = IndexR(i,j,l);
     1380       data_[ip] += snoise*corr*NorRand();
     1381     }
     1382   }
     1383 }
     1384 if(type_evol>0)
     1385   cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
     1386
     1387 if(intercor!=NULL) delete intercor;
     1388}
     1389
     1390}  // Fin namespace SOPHYA
     1391
     1392
     1393
     1394
     1395/*********************************************************************
    13311396void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
    13321397// Add AGN flux into simulation:
    … …  
    14231488 
    14241489}
    1425 
    1426 void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
    1427 // add noise to every pixels (meme les <=0 !)
    1428 // type_evol = 0 : pas d'evolution de la puissance du bruit
    1429 //             1 : evolution de la puissance du bruit avec la distance a l'observateur
    1430 //             2 : evolution de la puissance du bruit avec la distance du plan Z
    1431 //                 (tous les plans Z sont mis au meme redshift z de leur milieu)
    1432 {
    1433  if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
    1434  check_array_alloc();
    1435 
    1436  if(type_evol<0) type_evol = 0;
    1437  if(type_evol>2) {
    1438    char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
    1439    cout<<bla<<endl; throw ParmError(bla);
    1440  }
    1441 
    1442  vector<double> correction;
    1443  InterpFunc *intercor = NULL;
    1444 
    1445  if(type_evol>0) {
    1446    // Sigma_Noise(en mass) :
    1447    //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
    1448    //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
    1449    //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
    1450    //                                    nlobe ~ 1/Dtrcom^2
    1451    //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
    1452    if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
    1453      char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
    1454      cout<<bla<<endl; throw ParmError(bla);
    1455    }
    1456    InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
    1457    long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
    1458    for(long i=0;i<nsz;i++) {
    1459      double d = interpinv.X(i);
    1460      double zred = interpinv(d);
    1461      double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
    1462      double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
    1463      double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
    1464      double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
    1465      double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
    1466      if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
    1467        cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
    1468            <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
    1469      correction.push_back(corr);
    1470    }
    1471    intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
    1472  }
    1473 
    1474  double corrlim[2] = {1.,1.};
    1475  for(long i=0;i<Nx_;i++) {
    1476    double dx2 = DXcom(i); dx2 *= dx2;
    1477    for(long j=0;j<Ny_;j++) {
    1478      double dy2 = DYcom(j); dy2 *= dy2;
    1479      for(long l=0;l<Nz_;l++) {
    1480        double corr = 1.;
    1481        if(type_evol>0) {
    1482          double dz = DZcom(l);
    1483          if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
    1484            else dz = fabs(dz); // tous les plans Z au meme redshift
    1485          corr = (*intercor)(dz);
    1486          if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
    1487        }
    1488        int_8 ip = IndexR(i,j,l);
    1489        data_[ip] += snoise*corr*NorRand();
    1490      }
    1491    }
    1492  }
    1493  if(type_evol>0)
    1494    cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
    1495 
    1496  if(intercor!=NULL) delete intercor;
    1497 }
    1498 
    1499 }  // Fin namespace SOPHYA
     1490*********************************************************************/
Note: See TracChangeset for help on using the changeset viewer.