Changeset 4008 in Sophya

Timestamp:

Jun 23, 2011, 2:53:23 PM (13 years ago)

Author:

cmv

Message:

add a 2sd power law in frequency for GSM fit, cmv 23/06/2011

Location:

trunk/Cosmo/SimLSS

Files:

• cmvgsmcfit.cc (modified) (9 diffs)
• genefluct3d.cc (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvgsmcfit.cc

@@ -25 +25 @@
 public:
   MyFCN(TVector<r_4>& FreqMHz,TVector<r_4>& V,TVector<r_4>& EV);
+  void SetNu0(double nu0Mhz) {_nu0 = nu0Mhz;}
+  void SetTyp(bool law2 = false) {_law2 = law2;}
   void InitFit(MnUserParameters& upar);
   double Func(double fr,const vector<double>& par) const;
   double operator()(const vector<double>& par) const;
   double Up(void) const {return 1.;}
-  void SetNu0(double nu0Mhz) {_nu0 = nu0Mhz;}
   //protected:
   double _nu0;
+  bool _law2;
   mutable int _ndata;
   TVector<r_4> _F;
@@ -43 +45 @@
 {
   int nslicemax = -1;
+  int typfit = 0;
   if(narg<=1) {
-    cout<<"cmvgsmcfit cubefile.ppf [nslicemax]"<<endl;
+    cout<<"cmvgsmcfit cubefile.ppf [nslicemax] [typfit=0,1]"<<endl;
     return -1;
   }
 
   if(narg>2) nslicemax = atoi(arg[2]);
-  cout<<"Fitting "<<arg[1]<<" , nslicemax="<<nslicemax<<endl;
+  if(narg>3) typfit = atoi(arg[3]);
+  cout<<"Fitting "<<arg[1]<<" , nslicemax="<<nslicemax<<" , typfit="<<typfit<<endl;
 
 int rc = 0;
@@ -61 +65 @@
   pis>>PPFNameTag("freq")>>Freq;
 
-  const int nvar = 8;
-  const char *vname[nvar] = {"x2r","ok","ph","th","v0","vn","A","N"};
+  const int nvar = 10;
+  const char *vname[nvar] = {"x2r","ok","ph","th","v0","vn","A","N","B","DN"};
   NTuple nt(nvar,vname);
   double xnt[nvar]; for(int i=0;i<nvar;i++) xnt[i] = 0.;
@@ -70 +74 @@
     cout<<"...j = "<<j<<endl;
   for(int k=0;k<Cube.SizeZ();k++) {
+    //cout<<"...k = "<<k<<endl;
 
     TVector<r_4> V(Cube.SizeX()), EV(Cube.SizeX());
@@ -79 +84 @@
     MyFCN myfcn(Freq,V,EV);
     myfcn.SetNu0(Freq(0));
+    //myfcn.SetNu0((Freq(0)+Freq(Freq.Size()-1))/2.);
+    myfcn.SetTyp( ((typfit==1)? true: false) );
     MnUserParameters upar;
     myfcn.InitFit(upar);
@@ -102 +109 @@
     xnt[2] = j; xnt[3] = k;
     xnt[4] = V(0); xnt[5] = V(V.Size()-1);
-    for(unsigned int i=0;i<par.size();i++) xnt[6+i] = par[i];
+    for(int i=6;i<6+par.size() && i<nvar;i++) xnt[i] = par[i-6];
     nt.Fill(xnt);
     //vector<double> err = ustate.Errors();
@@ -155 +162 @@
     throw SzMismatchError("MyFCN::MyFCN_Error: incompatible F / V / EV sizes");
   _nu0 = _F(0);
+  _law2 = false;
   _ndata =_V.Size();
 }
@@ -165 +173 @@
 
   upar.Add("A",a,a/50.);
-  upar.Add("N",ni,fabs(ni/50.));
+  upar.Add("Nu",ni,fabs(ni/50.));
+
+  if(!_law2) return;
+  upar.Add("B",0.,a/200.);
+  upar.Add("DNu",0.,fabs(ni/200.));
+  
 }
 
 double MyFCN::Func(double fr,const vector<double>& par) const
-// v = A*(nu/nu0)^N
+// v = A*(nu/nu0)^Nu
+// v = A*(nu/nu0)^Nu + B/A*(nu/nu0)^(Nu+DNu) = A*(nu/nu0)^Nu * (1 + B*(nu/nu0)^DNu)
 {
  double f = par[0]*pow(fr/_nu0,par[1]);
+ if(_law2) f *= 1. + par[2]*pow(fr/_nu0,par[3]);
  return f;
 }
@@ -210 +225 @@
 n/plot nt.log10(v0/vn)%v0 ! ! "crossmarker3"
 
-
 zone
-n/plot nt.A%_nl ! ! "crossmarker3"
-n/plot nt.N%_nl ! ! "crossmarker3"
-
-n/plot nt.A%v0 ! ! "crossmarker3"
-n/plot nt.(A-v0)/v0%_nl ! ! "crossmarker3"
-n/plot nt.N%A ! ! "crossmarker3"
+n/plot nt.A%_nl ok>0 ! "crossmarker3"
+n/plot nt.N%_nl ok>0 ! "crossmarker3"
+n/plot nt.B%_nl ok>0 ! "crossmarker3"
+n/plot nt.DN%_nl ok>0 ! "crossmarker3"
+
+n/plot nt.A%v0 ok>0 ! "crossmarker3"
+n/plot nt.(A-v0)/v0%_nl ok>0 ! "crossmarker3"
+
+n/plot nt.N%A ok>0 ! "crossmarker3"
+n/plot nt.B%A ok>0 ! "crossmarker3"
+n/plot nt.DN%N ok>0 ! "crossmarker3"
 
 set n 0

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -964 +964 @@
 ---- ATTENTION: Le spectre Pk doit etre (dRho/Rho)(k)
 ---- On genere la vitesse particuliere (dans le referentiel comobile au redshift "z")
-     Vlos(k) = i*Beta(z)*k(los)/k^2*(dRho/Rho)(k)
+     Vlos(k) = i*Beta(z)*k(los)/k^2*(dRho/Rho)(k)    (unite Mpc)
         .avec Beta(z) = dln(D)/da = -(1+z)/D*dD/dz
-        .on convertit en vitesse en faisant Vlos(k)/H(z)
+        .puis on convertit en vitesse en faisant Vlos(k) * H(z)    (unite lm/s)
 */
 {