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

#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;
 int_8 ntfill = 100000;
 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 naxis1 = f3dr.ReadKeyL("NAXIS1");
 long naxis2 = f3dr.ReadKeyL("NAXIS2");
 long naxis3 = f3dr.ReadKeyL("NAXIS3");
 cout<<"Naxis: 1="<<naxis1<<" 2="<<naxis2<<" 3="<<naxis3<<endl;
 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);

 const int nvar = 3;
 const char *vname[nvar] = {"l","ll","d"};
 NTuple nt(nvar,vname);
 r_4 xnt[nvar];
 int_8 ntmod = fluct3d.NPix()/ntfill; if(ntmod==0) ntmod = 1;

 cout<<"> filling redshift distorted cube"<<endl;
 int_8 nf = 0;
 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);
     int ll = int((double)l + d/dz + 0.5);
     if(ll<0 || ll>=nz) continue;
     Rdis(ll) += R(l);
     if(ntfill>0 && nf%ntmod==0)
       {xnt[0]=l; xnt[1]=ll; xnt[2]=d; nt.Fill(xnt);}
     nf++;
   }
   for(int l=0;l<nz;l++) rgen(l,j,i) += Rdis(l);
 }
 }
 PrtTim(">>>> End filling redshift distorted cube");
 pos.PutObject(nt,"nt");

 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

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

 */