source: Sophya/trunk/Cosmo/SimLSS/cmvrvloscor.cc@ 3771

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "sophyainit.h"
#include "timing.h"
#include "dvlist.h"
#include "ntuple.h"
#include "fabtcolread.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "genefluct3d.h"

void usage(void);
void usage(void)
{
  cout<<"cmvrvloscor rho.fits vlos.fits"<<endl;
}

int main(int narg,char *arg[])
{
 int nthread = 1;
 SophyaInit();
 InitTim();

 if(narg<=2) {usage(); return -1;}

 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
 try {
 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH

 cout<<"> read rho: "<<arg[1]<<endl;
 FitsImg3DRead f3dr(arg[1],0,5);
 cout<<"> read vlos: "<<arg[2]<<endl;
 FitsImg3DRead f3dv(arg[2],0,5);
 long Nx = f3dr.ReadKeyL("Nx");
 long Ny = f3dr.ReadKeyL("Ny");;
 long Nz = f3dr.ReadKeyL("Nz");;
 cout<<"N: x="<<Nx<<" y="<<Ny<<" z="<<Nz<<endl;
 double Dx = f3dr.ReadKey("Dx");
 double Dy = f3dr.ReadKey("Dy");
 double Dz = f3dr.ReadKey("Dz");
 cout<<"D: x="<<Dx<<" y="<<Dy<<" z="<<Dz<<endl;
 double Zref = f3dr.ReadKey("ZREF");
 double Href = f3dr.ReadKey("HREF");
 cout<<"Zref="<<Zref<<" Href="<<Href<<endl;

 POutPersist pos("cmvrvloscor.ppf");

 GeneFluct3D fluct3d(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
 fluct3d.Print();
 TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
 rgen = 0.;
 TVector<r_8> R(Nz), Rdis(Nz);

 cout<<"> filling redshift distorted cube"<<endl;
 for(int i=0;i<Nx;i++) {
   if(i%(Nx/10)==0) cout<<"i="<<i<<endl;
 for(int j=0;j<Ny;j++) {
   for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
   Rdis = 0.;
   for(int l=0;l<Nz;l++) {
     double d = (1.+Zref) / Href * f3dv.Read(l,j,i);
     double lpd = (double)l + d/Dz; // valeur du deplacee
     // on repartit proportionnelement au recouvrement sur 2 pixels
     int l1 = int(lpd); // pixel de droite
     int l2 = l1 + 1;  // pixel de gauche
     lpd -= (double)l1;  // recouvrement du pixel du dessus
     if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
     if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
   }
   for(int l=0;l<Nz;l++) rgen(l,j,i) += Rdis(l);
 }
 }
 PrtTim(">>>> End filling redshift distorted cube");

 fluct3d.ReComputeFourier();
 PrtTim(">>>> End ReComputing spectrum");

 cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
 double dkmin = fluct3d.GetKincMin();
 double knyqmax = fluct3d.GetKmax();
 long nherr = long(knyqmax/dkmin+0.5);
 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
 HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
 hpkrec.ReCenterBin(); hpkrec.Show();
 fluct3d.ComputeSpectrum(hpkrec);
 pos.PutObject(hpkrec,"hpkrec");
 PrtTim(">>>> End Computing final spectrum");

 cout<<"\n--- Computing final 2D spectrum"<<endl;
 double dktmin = fluct3d.GetKTincMin();
 double ktnyqmax = fluct3d.GetKTmax();
 long nherrt = long(ktnyqmax/dktmin+0.5);
 double dkzmin = fluct3d.GetKinc()[2];
 double kznyqmax = fluct3d.GetKnyq()[2];
 long nherrz = long(kznyqmax/dkzmin+0.5);
 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
     <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
 hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
 fluct3d.ComputeSpectrum2D(hpkrec2);
 pos.PutObject(hpkrec2,"hpkrec2");
 PrtTim(">>>> End Computing final 2D spectrum");

 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
 } catch (PException& exc) {
   cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
   return 77;
 } catch (std::exception& sex) {
   cerr << "cmvrvloscor.cc std::exception :" 
        << (string)typeid(sex).name() << "\n msg= " 
        << sex.what() << endl;
   return 78;
 } catch (...) {
   cerr << "cmvrvloscor.cc catched unknown (...) exception  " << endl; 
   return 79;
 }
 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH

 return 0;
}

/*
openppf cmvrvloscor.ppf

n/plot hpkrec.val%x x>0 ! "nsta cpts logx"

imag hpkrec2
addoval 0 0 0.05 0.05 "green" false
addoval 0 0 0.1 0.1 "green" false
addoval 0 0 0.25 0.25 "green" false
addoval 0 0 0.5 0.5 "green" false
x = ${hpkrec2.xmax} / 2.
addoval 0 0 $x $x "green" false
x = ${hpkrec2.ymax} / 2.
addoval 0 0 $x $x "green" false

# proj selon kT (black), selon kZ (red)
n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta cpts"
 */