Changeset 3313 in Sophya for trunk/Cosmo

Timestamp:

Aug 24, 2007, 4:59:12 PM (18 years ago)

Author:

cmv

Message:

nettoyage et mise en forme cmv 24/08/2007

File:

• trunk/Cosmo/SimLSS/cmvtuniv.cc (modified) (13 diffs)

trunk/Cosmo/SimLSS/cmvtuniv.cc

@@ -13 +13 @@
 #include "geneutils.h"
 
-void usage(void);
-void usage(void) {cout<<"cmvtuniv z1,z2,dz [perc,dzinc,dzmax,glorder]"<<endl;}
-
 void tstprint(CosmoCalc& univ,double z1,double z2,double dz);
 void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz);
@@ -23 +20 @@
 
 const double deg2rad = M_PI/180.;
+
+void usage(void);
+void usage(void) {
+ cout
+ <<"cmvtuniv [options] z1,z2"<<endl
+ <<"  -D dz"<<endl
+ <<"  -U h100,om0,or0,ol0,w0"<<endl
+ <<"  -F flat(0,1,2)"<<endl
+ <<"  -I perc,dzinc,dzmax,glorder"<<endl
+ <<" --- What has to be done:"<<endl
+ <<"  -W 1/2/3/4/5"<<endl
+ <<"     1 : print info from z1 to z2 by step dz"<<endl
+ <<"     2 : compare with/without spline from z1 to z2 by step dz"<<endl
+ <<"     3 : test speed from z1 to z2 by step dz"<<endl
+ <<"     4 : fill ntuple with infos from z1 to z2 by step dz"<<endl
+ <<"     5 : test interpolation loscom from z1 to z2 by step dz"<<endl
+ <<endl;
+}
 
 int main(int narg,char *arg[])
@@ -36 +51 @@
  //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.; flat = 2;
 
+ // -- parametre d'integration
+ double perc=0.01,dzinc=-1.,dzmax=-1.;
+ unsigned short glorder=4;
+ // -- redshift
  double z1=0., z2=-1., dz=-1.;
- if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);
+ // -- Work to be done
+ int todo = 1;
+
+ InitTim();
+
+ // Read arguments
+ char c;
+ while ((c = getopt(narg, arg, "hD:U:F:I:W:")) != -1) {
+  switch (c) {
+  case 'D' :
+    sscanf(optarg,"%lf",&dz);
+    break;
+  case 'U' :
+    sscanf(optarg,"%lf,%lf,%lf,%lf,%lf",&h100,&om0,&or0,&ol0,&w0);
+    break;
+  case 'F' :
+    sscanf(optarg,"%hu",&flat);
+    break;
+  case 'I' :
+    sscanf(optarg,"%lf,%lf,%lf,%hu",&perc,&dzinc,&dzmax,&glorder);
+    break;
+  case 'W' :
+    sscanf(optarg,"%d",&todo);
+    break;
+  case 'h' :
+  default:
+    usage();
+    return -1;
+  }
+}
+ if(optind<narg) sscanf(arg[optind],"%lf,%lf",&z1,&z2);
+
+ if(z1<0.) z1 = 0.;
  if(z2<=z1) z2 = z1+1.;
- if(dz<0.) dz = 10.*(z2-z1); 
+ if(dz<=0.) dz = 10.*(z2-z1);
  cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
-
- double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
- if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%hu",&perc,&dzinc,&dzmax,&glorder);
- cout<<" perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
-
- CosmoCalc univ(flat,true,z2);
- univ.SetInteg(perc,dzinc,dzmax,glorder);
- univ.SetDynParam(h100,om0,or0,ol0,w0);
- univ.PrtInteg();
- univ.Print();
-
- CosmoCalc univ1(univ);
+ cout<<"perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
+ cout<<"h100="<<h100<<" om0="<<om0<<" or0="<<or0<<" ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
+ cout<<"todo ="<<todo<<endl;
+
+ CosmoCalc univ1(flat,true,z2);
  univ1.SetInteg(perc,dzinc,dzmax,glorder);
  univ1.SetDynParam(h100,om0,or0,ol0,w0);
  univ1.PrtInteg();
  univ1.Print();
-
- tstprint(univ,z1,z2,dz);
- //tstprint(univ,univ1,z1,z2,dz);
- //tstspeed(univ,z1,z2,dz);
- //tstntuple(univ,z1,z2,dz);
- //tstinterp(univ,z1,z2,dz);
+ cout<<endl;
+
+ if(todo==1)  tstprint(univ1,z1,z2,dz);
+ else if(todo==2) {
+   CosmoCalc univ2(flat,false,z2);
+   univ2.SetInteg(perc,dzinc,dzmax,glorder);
+   univ2.SetDynParam(h100,om0,or0,ol0,w0);
+   univ2.PrtInteg();
+   univ2.Print();
+   tstprint(univ1,univ2,z1,z2,dz);
+ }
+ else if(todo==3) tstspeed(univ1,z1,z2,dz);
+ else if(todo==4) tstntuple(univ1,z1,z2,dz);
+ else if(todo==5)tstinterp(univ1,z1,z2,dz);
+ else cout<<"Unknown job to be done !!!!!  todo="<<todo<<endl;
 
  return 0;
@@ -70 +123 @@
 void tstprint(CosmoCalc& univ,double z1,double z2,double dz)
 {
- cout<<"\nTSTPRINT()"<<endl;
- double dzz = dz; if(dzz>z2-z1) dzz = 0.01;
- for(double z=z1;z<z2+dz/2.;z+=dz) {
-   cout<<"--"<<endl;
+ if(dz>z2-z1) dz = z2-z1;
+ cout<<"\nTSTPRINT(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   cout<<"-- z = "<<z<<endl;
    univ.Print(z);
-   cout<<"Volume comoving in [z,z+"<<dzz<<"] for 4Pi sr: "
-       <<univ.Vol4Pi(z,z+dzz)<<" Mpc^3"<<endl;
+   cout<<"Volume comoving in [z,z+"<<dz<<"] for 4Pi sr: "
+       <<univ.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
    double dang = univ.Dang(z);
    double a = deg2rad/3600.;
@@ -86 +139 @@
    double dloscom = univ.Dhubble() / univ.E(z);
    cout<<"dz=1 -> "<<dloscom<<" Mpc com"<<endl;
-
- }
+ }
+
 }
 
@@ -93 +146 @@
 void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz)
 {
- cout<<"\nTSTPRINT()"<<endl;
- for(double z=z1;z<z2+dz/2.;z+=dz) {
-   cout<<endl<<"--spline:"<<endl;
+ if(dz>z2-z1) dz = z2-z1;
+ cout<<"\nTSTPRINT(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   cout<<endl<<"--spline: z = "<<z<<endl;
    univ1.Print(z);
    cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
        <<univ1.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
-   cout<<"--integ:"<<endl;
+   cout<<"--integ: z = "<<z<<endl;
    univ2.Print(z);
    cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
@@ -109 +163 @@
 void tstspeed(CosmoCalc& univ,double z1,double z2,double dz)
 {
- cout<<"\nTSTSPEED()"<<endl;
+ if(dz>z2-z1) dz = z2-z1;
+ cout<<"\nTSTSPEED(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
  double sum = 0.;
  int_4 n=0;
+ PrtTim("-- Beginning");
  for(double z=z1;z<z2+dz/2.;z+=dz) {
    sum += univ.Dang(z);
    n++;
  }
- if(n>0) cout<<"... "<<sum/n<<endl;
+ if(n>0) cout<<"n="<<n<<" ...dum="<<sum/n<<endl;
+ PrtTim("-- Ending");
 }
 
@@ -122 +179 @@
 void tstntuple(CosmoCalc& univ,double z1,double z2,double dz)
 {
- cout<<"\nTSTNTUPLE()"<<endl;
+ if(dz>z2-z1) dz = z2-z1;
+ cout<<"\nTSTNTUPLE(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
  double norm =  univ.Dhubble();  // 1.
  double norm3 = pow(norm,3.);
@@ -152 +210 @@
  }
  cout<<">>>> Ecriture"<<endl;
- string tag = "cmvtuniv.ppf";
+ string tag = "cmvtuniv_nt.ppf";
  POutPersist pos(tag);
  tag = "nt"; pos.PutObject(nt,tag);
@@ -159 +217 @@
 
 /*
-openppf cmvtuniv.ppf
+openppf cmvtuniv_nt.ppf
 
 set cut 1
@@ -191 +249 @@
 void tstinterp(CosmoCalc& univ,double z1,double z2,double dz)
 {
- cout<<"\nTSTINTERP()"<<endl;
-
- // On repmplit les donnes completes
+ if(dz>z2-z1) dz = z2-z1;
+ cout<<"\nTSTINTERP(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
+ // On remplit les donnes completes
  vector<double> Z,D;
  for(double z=z1;z<z2+dz/2.;z+=dz) {
@@ -226 +284 @@
  }
  cout<<">>>> Ecriture"<<endl;
- string tag = "cmvtuniv.ppf";
+ string tag = "cmvtuniv_int.ppf";
  POutPersist pos(tag);
  tag = "nt"; pos.PutObject(nt,tag);
@@ -232 +290 @@
 
 /*
-openppf cmvtuniv.ppf
+openppf cmvtuniv_int.ppf
 
 n/plot nt.d%z ! ! "nsta connectpoints"