Changeset 3932 in Sophya for trunk/Cosmo

Timestamp:

Dec 23, 2010, 7:33:37 PM (15 years ago)

Author:

ansari

Message:

amelioration repicon.cc (decodage arguments) , Reza 23/12/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• mdish.cc (modified) (2 diffs)
• mdish.h (modified) (2 diffs)
• repicon.cc (modified) (6 diffs)

trunk/Cosmo/RadioBeam/mdish.cc

@@ -183 +183 @@
   SetThetaPhiRange();
   SetRespHisNBins();
+  SetPrtLevel();
   mcnt_=0;
 }
@@ -207 +208 @@
 
   double dtet = thetamax_/(double)ntet_;
-  double dphi = phimax_/(double)ntet_;
+  double dphi = phimax_/(double)nphi_;
+  cout << " MultiDish::GetResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_ 
+       << " PhiMax=" <<  phimax_ << " NPhi=" << nphi_ << endl;
 
   double sumw = 0.;
-  for(int kt=0; kt<ntet_; kt++) 
+  for(int kt=0; kt<ntet_; kt++) {
+    double theta=(double)kt*dtet;
     for(int jp=0; jp<nphi_; jp++) {
-      sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi);
-      sumw += CumulResp(rd, (double)kt*dtet, -(double)jp*dphi);
-      sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi+M_PI);
-      sumw += CumulResp(rd, (double)kt*dtet, -(double)jp*dphi-M_PI);
+      double phi=(double)jp*dphi;
+      sumw += CumulResp(rd, theta, phi);
+      if (theta<1.e-9) continue;
+      sumw += CumulResp(rd, theta, -phi);
+      sumw += CumulResp(rd, theta, phi+M_PI);
+      sumw += CumulResp(rd, theta, -(phi+M_PI));
     }
+    if (prtlev_>0) cout << "  MultiDish::GetResponse() done ktheta=" << kt << " / MaxNTheta= " << ntet_ << endl; 
+  }
   double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
   double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;

trunk/Cosmo/RadioBeam/mdish.h

@@ -120 +120 @@
   MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto=false);
 
+  // Pour phi, les angles phi, -phi, phi+pi, -(phi+pi) sont prises en compte
   inline void SetThetaPhiRange(double thetamax=0., int ntet=1, double phimax=0., int nphi=1)
   { thetamax_=thetamax; ntet_=ntet; phimax_=phimax; nphi_=nphi; }
+
+  inline int SetPrtLevel(int lev=0, int prtmod=10) 
+  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
 
   inline void SetRespHisNBins(int nx=128, int ny=128)
@@ -145 +149 @@
   QHis2D h2w_;
   int mcnt_;
+  int prtlev_,prtmodulo_;
 };
 

trunk/Cosmo/RadioBeam/repicon.cc

@@ -3 +3 @@
     R. Ansari , C. Magneville - Juin-Dec 2010 
 
-  Usage:  repicon configId OutPPFName [z_Redshift=0.7] [RenormalizeMax]
+  Usage:  repicon [-parname value] configId OutPPFName 
 ---------------------------------------------------------------  */
 
@@ -39 +39 @@
 // ---------------------------------------------------------------------
 
+void Usage() 
+{
+  cout << " Usage: repicon [-parname Value] configId OutPPFName \n" 
+       << " configIds: f4x4,f8x8,f20x20, confA,confB,confC,confD, hex12,cross11, f3cyl,f8cyl,f3cylp,f8cylp \n"
+       << "   f4x4 , f8x8 , f20x20 Filled array of nxn dishes  \n" 
+       << "   confA , confB, confC, confD : semi-filled array of dishes \n"
+       << "   hex12,cross11 : ASKAP like double hexagonal (12xD=12m), cross config (11xD=12m) \n"
+       << "   f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders \n"
+       << "  [ -parname value] : -renmax -z -prt -D -lmax \n" 
+       << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
+       << "    -z redshift (default=0.7) --> determines Lambda \n"
+       << "    -D DishDiameter (default=5 m) \n"
+       << "    -lmax array extension (default=100 m ) for response calculation kmax \n" 
+       << "    -rot ThetaMaxDeg,NTheta,PhiMaxDeg,NPhi (default NO-rotate/pointing -> 23,10,30,15 ) \n" 
+       << "    -prt PrtLev,PrtModulo (default=0,10) \n" 
+       << endl; 
+  return;
+}
+
 //      ------------------ MAIN PROGRAM ------------------------------
 int main(int narg, const char* arg[])
 {
   if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
-    cout << " Usage: repicon configId OutPPFName [z_redshift=0.7] [RenormalizeMax] \n" 
-         << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes (D=5m) \n" 
-         << "          confA , confB, confC, confD : semi-filled array of dishes \n"
-         << "          hex12,cross11 : ASKAP like double hexagonal (12xD=12m), cross config (11xD=12m) \n"
-         << "          f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders " << endl; 
+    Usage();
     return 1;
   }
-  cout << " ==== repicon.cc program , (u,v) plane response ==== " << endl;
   // make sure SOPHYA modules are initialized 
   SophyaInit();  
@@ -57 +71 @@
   //--- decoding command line arguments 
   string config = "f8x8" ;
-  if (narg>1) config = arg[1];
   string outfile = "repicon.ppf";  
-  if (narg>2) outfile = arg[2];
-  if (outfile==".")  outfile = "repicon.ppf";
-  double LAMBDA=0.357 ;  // 21 cm at z=0.7 
-  if (narg>3) LAMBDA = atof(arg[3]);
   bool fgrenorm=false;
   double rmax=1.;
-  if (narg>4) {
-    rmax=atof(arg[4]);
-    fgrenorm=true;
-  }
+  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
+  int prtlev=0;
+  int prtmod=10;
+  
+  double Ddish=5.;
+  bool fgDfixed=false;
+  double LMAX = 100.;  // taille de la zone, pour calcul kmax 
+  bool fgLMAXfixed=false;
+  double thetamxdeg=23.;  // 23 degres : angle d'inclinaison de l'orbite terrestre
+  double phimxdeg=30.;
+  int nteta=10;
+  int nphi=15;
+  bool fgpoint=false;
+
+  int ka=1;
+  while (ka<(narg-1)) {
+    if (strcmp(arg[ka],"-renmax")==0) {
+      rmax=atof(arg[ka+1]);  fgrenorm=true;   ka+=2;
+    }
+    else if (strcmp(arg[ka],"-z")==0) {
+      z_Redshift=atof(arg[ka+1]);  ka+=2;
+    }
+    else if (strcmp(arg[ka],"-D")==0) {
+      Ddish=atof(arg[ka+1]);  fgDfixed=true;  ka+=2;
+    }
+    else if (strcmp(arg[ka],"-lmax")==0) {
+      LMAX=atof(arg[ka+1]);  fgLMAXfixed=true;  ka+=2;
+    }
+    else if (strcmp(arg[ka],"-rot")==0) {
+      sscanf(arg[ka+1],"%lg,%d,%lg,%d",&thetamxdeg,&nteta,&phimxdeg,&nphi);  fgpoint=true;  ka+=2;
+    }
+    else if (strcmp(arg[ka],"-prt")==0) {
+      sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod);   ka+=2;
+    }
+    else break; 
+  }
+
+  if ((ka+1)>=narg) {
+    cout << " repicon / Argument error " << endl;
+    Usage();
+    return 2;
+  }
+
+  config = arg[ka];
+  outfile = arg[ka+1];
   //-- end command line arguments
   
   int rc = 1;  
   try {  // exception handling try bloc at top level
-
-    double Ddish = 5.;
-    double Ddish2 = 7.5;
+    cout << " ==== repicon.cc program , (u,v) plane response ==== " << endl;
+    
     double Eta=0.95;
     double cylW=12.;   // Largeur des cylindres 
@@ -82 +131 @@
     double etaRL=0.9;   // Efficacite de couverture le long du cylindre
 
-    double D = 100.;    // Taille de la zone
-
     int cnt=0;
 
@@ -89 +136 @@
     cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config << endl;
 
-    bool fgpoint=false;
     if (config=="f4x4") {
       vdishes=CreateFilledSqConfig(4,Ddish, Eta);
@@ -116 +162 @@
     }
     else if (config=="confA") {
-      fgpoint=true;
       vdishes=CreateConfigA(Ddish, Eta);
     }
     else if (config=="confB") {
-      fgpoint=true;
       vdishes=CreateConfigB(Ddish, Eta);
     }
     else if (config=="confC") {
-      fgpoint=true;
       vdishes=CreateConfigC(Ddish, Eta);
     }
     else if (config=="confD") {
-      fgpoint=true;
       vdishes=CreateConfigD(Ddish, Eta);
     }
     else if (config=="cross11") {
-      fgpoint=true;
-      Ddish = 12.;
+      if (!fgDfixed) Ddish = 12.;
       double base=20.;
       Eta=0.95;
-      D=250.;
+      if (!fgLMAXfixed) LMAX = 250.;
       vdishes=CreateCrossConfig(Ddish,base,Eta);
     }
     else if (config=="hex12") {
-      fgpoint=true;
-      Ddish = 12.;
+      if (!fgDfixed) Ddish = 12.;
       Eta=0.95;
-      D=350.;
+      if (!fgLMAXfixed) LMAX = 350.;
       vdishes=CreateDoubleHexagonConfig(Ddish);
     }
-
     else {
       cout << " NON valid configuration option -> exit" << endl;
       return 99;
     }
-
-    double Dol = D/LAMBDA;
-    double LMAX = D;
+    
+    H21Conversions conv;
+    double LAMBDA=0.357 ;  // 21 cm at z=0.7 
+    conv.setRedshift(z_Redshift);
+    LAMBDA = conv.getLambda();
+    
+    double Dol = LMAX/LAMBDA;
     bool fgnoauto = true;
     int NRX=200;
     int NRY=200;
-
+    
+    cout << " repicon[1] : Lambda=" << LAMBDA << " LMAX= " << LMAX << " NDishes=" << vdishes.size() 
+         << " D-Dish=" << Ddish << " m." << endl;
+ 
     MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
+    mdish.SetPrtLevel(prtlev,prtmod);
     mdish.SetRespHisNBins(NRX,NRY); 
-    if (fgpoint)  // 23 degres : angle d'inclinaison de l'orbite terrestre
-      mdish.SetThetaPhiRange(Angle(23.,Angle::Degree).ToRadian(),10, M_PI/6., 15);
+
+    if (fgpoint)  {
+      cout << " repicon[1.b] : activating pointing , ThetaMaxDeg=" << thetamxdeg << " NTheta=" << nteta
+           << "  PhiMaxDeg= " << phimxdeg << " NPhi" << nphi << endl;
+      mdish.SetThetaPhiRange(Angle(thetamxdeg,Angle::Degree).ToRadian(),nteta,Angle(phimxdeg,Angle::Degree).ToRadian(), nphi);
+    }
     cout << " repicon[2] : calling mdish.GetResponse() ..."<< endl;