Changeset 3199 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.cc

Timestamp:

Apr 5, 2007, 7:19:06 PM (18 years ago)

Author:

cmv

Message:

Intro raffinement dans simul AGN, variation dNu Da Dlum selon Oz cmv 05/04/2007

File:

• trunk/Cosmo/SimLSS/cmvobserv3d.cc (modified) (14 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -31 +31 @@
      <<" -2 : compute 2D spectrum"<<endl
      <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
-     <<" -A <log10(S_agn en Msol)>,sigma : distribution des AGN par pixel"<<endl
+     <<" -A <log10(S_agn en Jy)>,sigma,powlaw : distribution des AGN par pixel"<<endl
+     <<" -K : on utilise le calcul de spectre bricole pour les AGN (bidon!)"<<endl
      <<" -W : write cube in FITS format"<<endl
      <<" -P : write cube in PPF format"<<endl
@@ -79 +80 @@
  // *** AGN
  bool do_agn = false;
- double lmsol_agn=-99., lsigma_agn=0.;   // en equivalent MSol
+ double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
+ double powlaw_agn = 0.;
+ bool killkz = false;
 
  // *** type de generation
@@ -95 +98 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2GKx:y:z:s:Z:M:A:")) != -1) {
   switch (c) {
   case 'a' :
@@ -129 +132 @@
   case 'A' :
     do_agn = true;
-    sscanf(optarg,"%lf,%lf",&lmsol_agn,&lsigma_agn);
+    sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
+    break;
+  case 'K' :
+    killkz = true;
     break;
   case 'V' :
@@ -161 +167 @@
      <<", schnpt="<<schnpt<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
- if(do_agn) cout<<"AGN: <log10(Msol)>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
+ if(do_agn)
+   cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+       <<" , powlaw="<<powlaw_agn<<endl;
 
  //-----------------------------------------------------------------
@@ -271 +279 @@
  fluct3d.SetCosmology(univ);
  fluct3d.SetGrowthFactor(growth);
- fluct3d.LosComRedshift(0.001);
+ fluct3d.LosComRedshift(0.001,-1);
  TArray<r_8>& rgen = fluct3d.GetRealArray();
  cout<<endl; fluct3d.Print();
@@ -383 +391 @@
    cout<<"\n--- Apply Growth factor"<<endl;
    cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
-   fluct3d.ApplyGrowthFactor(-1);
+   fluct3d.ApplyGrowthFactor();
    rmin,rmax; rgen.MinMax(rmin,rmax);
    cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
@@ -486 +494 @@
 
  if(do_agn) {
-   cout<<"\n--- Add AGN: <mass>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
-   fluct3d.AddAGN(lmsol_agn,lsigma_agn);
+   cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+       <<" , powlaw="<<powlaw_agn<<endl;
+   fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
    nm = fluct3d.MeanSigma2(rm,rs2);
    cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
@@ -503 +512 @@
  }
 
+ if(wfits) {
+   fluct3d.WriteFits("!cmvobserv3d_rf.fits");
+   PrtTim(">>>> End WriteFits");
+ }
+ if(wppf) {
+   fluct3d.WritePPF("cmvobserv3d_rf.ppf",true);
+   PrtTim(">>>> End WritePPF");
+ }
+
  //-----------------------------------------------------------------
  // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
@@ -526 +544 @@
    hpkrec.ReCenterBin();
    hpkrec.Show();
-   fluct3d.ComputeSpectrum(hpkrec);
+   if(killkz) fluct3d.ComputeSpectrum_bricolo(hpkrec);
+     else     fluct3d.ComputeSpectrum(hpkrec);
    tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
    PrtTim(">>>> End Computing final spectrum");
@@ -536 +555 @@
    hpkrec2.ReCenterBin(); hpkrec2.Zero();
    hpkrec2.Show();
-   fluct3d.ComputeSpectrum2D(hpkrec2);
+   if(killkz) fluct3d.ComputeSpectrum2D_bricolo(hpkrec2);
+     else     fluct3d.ComputeSpectrum2D(hpkrec2);
    {
    tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
@@ -553 +573 @@
 readfits cmvobserv3d_r0.fits
 readfits cmvobserv3d_r.fits
+readfits cmvobserv3d_rf.fits
 
 openppf cmvobserv3d_k0.ppf
@@ -558 +579 @@
 openppf cmvobserv3d_r0.ppf
 openppf cmvobserv3d_r.ppf
-
-# pour le plot 2D d'une slice en Z du 3D: to2d nom_obj3D num_slice
-defscript to2d
+openppf cmvobserv3d_rf.ppf
+
+# pour le plot 2D d'une slice en Z du 3D: xy2d nom_obj3D num_slice
+defscript xy2d
  objaoper $1 sliceyz $2
- mv sliceyz_${2} ${1}_$2
- disp ${1}_$2
- echo display slice $2 of $1
+ mv sliceyz_${2} ${1}_Z_$2
+ disp ${1}_Z_$2
+ echo display slice $2 of $1 name is ${1}_Z_$2
+endscript
+ 
+# pour le plot 2D d'une slice en Y du 3D: xz2d nom_obj3D num_slice
+defscript xz2d
+ objaoper $1 slicexy $2
+ mv slicexy_${2} ${1}_Y_$2
+ disp ${1}_Y_$2
+ echo display slice $2 of $1 name is ${1}_Y_$2
 endscript 
-
-to2d $cobj 0
+ 
+# pour le plot 2D d'une slice en X du 3D: yz2d nom_obj3D num_slice
+defscript yz2d
+ objaoper $1 slicexz $2
+ mv slicexz_${2} ${1}_X_$2
+ disp ${1}_X_$2
+ echo display slice $2 of $1 name is ${1}_X_$2
+endscript 
+
+xy2d $cobj 0
+xz2d $cobj 0
+yz2d $cobj 0
 
 ######################################################
@@ -575 +615 @@
 set k pow(10.,x)
 n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
+
+echo ${hpkgen.sum}
+echo ${hpkgenf.sum}
+echo ${hpkrec.sum}
 
 zone