Changeset 2051 in Sophya


Ignore:
Timestamp:
Jun 11, 2002, 10:39:57 AM (23 years ago)
Author:
ansari
Message:

Ajout save des Alm ds map2cl.cc - Reza 11/6/2002

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaProg/PrgMap/map2cl.cc

    r1841 r2051  
    3434   -iter_order lval : 1,2,3,4... order of an iterative analysis,
    3535              3rd order is usually optimal (default=0 -> standard analysis)
    36    -fitsout: Select the FITS format for the output map (default PPF format)
    37    -fitsin : Select the FITS format for the input vector (default PPF format)
     36   -keepalm : Save Alm's as a 2-D array (Matrix)
     37   -fitsout: Select the FITS format for the output vector/alm's (default PPF format)
     38   -fitsin : Select the FITS format for the input map (default PPF format)
    3839   InFileName : Input file name (HEALPix map)
    39    OutFileName : Output file name (the C_l vector)
     40   OutFileName : Output file name (the C_l vector / and Alm's)
    4041
    4142  \endverbatim
     
    5354    cout << " map2cl : Spherical harmonics analysis - HEALPix map -> Power spectrum C_l \n"
    5455         << " Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg] [-iter_order lval]\n"
    55          << "        [-fitsin] [-fitsout] InFileName OutFileName \n"
     56         << "        [-keepalm AlmFileName] [-fitsin] [-fitsout] InFileName OutFileName \n"
    5657         << "   -float (-r4): single precision C_l and map (default = double) \n"
    5758         << "   -lmax lval: Maximum value for l index (default=255)\n"
    5859         << "   -thetacut dtdeg : Symmetric delta-theta cut (in degree) along equator \n"
    5960         << "                    (default=0 -> no cut)\n"
    60          << "   -iter_order lval : 1,2,3,4... order of an iterative analysis , 3rd order is usually optimal (default=0 -> standard analysis)\n"
    61          << "   -fitsout: Select the FITS format for the output map (default PPF format) \n"
    62          << "   -fitsin : Select the FITS format for the input vector (default PPF format) \n"
     61         << "   -iter_order lval : 1,2,3,4... order of an iterative analysis\n"
     62         << "       3rd order is usually optimal (default=0 -> standard analysis)\n"
     63         << "   -keepalm : Save Alm's as a 2-D complex array (Matrix<complex>) - only with PPF format  \n"
     64         << "   -fitsout: Select the FITS format for the output files (Map/Alm) (default PPF format) \n"
     65         << "   -fitsin : Select the FITS format for the input map (default PPF format) \n"
    6366         << "   InFileName : Input file name (HEALPix map) \n"
    6467         << "   OutFileName : Output file name (the C_l vector) \n" << endl;
     
    7275public :
    7376static void ComputeCl(string & infile, string & outfile, int lmax, double tcut,
    74                       int iterationOrder, bool fgfitsin, bool fgfitsout)
     77                      int iterationOrder, bool fgfitsin, bool fgfitsout, bool keepalm)
    7578{
    7679  double deg2rad =  M_PI/180.;
     
    98101  // Decomposition de la carte en C_l
    99102  SphericalTransformServer<T> sphtr;
    100   TVector<T> clvec = sphtr.DecomposeToCl(sph, lmax, ctcut, iterationOrder);
    101  
     103  TVector<T> clvec;
     104  TMatrix< complex<T> > almmtx;
     105  if (keepalm) {
     106    Alm<T> alm(lmax);
     107    sphtr.DecomposeToAlm(sph, alm, lmax, ctcut, iterationOrder);
     108    almmtx.SetSize(alm.Lmax()+1, alm.Lmax()+1);
     109    for(int il=0; il<alm.Lmax()+1; il++)
     110      for(int im=0; im<il; im++)
     111        almmtx(il, im) = alm(il, im);
     112    clvec = alm.powerSpectrum();
     113  }
     114  else {
     115    clvec = sphtr.DecomposeToCl(sph, lmax, ctcut, iterationOrder);
     116  }
    102117  T min, max;
    103118  double mean, sigma;
     
    109124 
    110125  if (fgfitsout) {
     126    FitsOutFile fio(outfile, FitsFile::clear);
    111127    cout << "--- Writing C_l vector to Output FITS file " << outfile << endl;
    112     FitsOutFile fio(outfile, FitsFile::clear);
    113128    fio << clvec;
    114129  }
    115130  else {
    116     cout << "--- Writing C_l vector to Output PPF file " << outfile << endl;
    117131    POutPersist ppo(outfile);
    118     ppo << clvec;
     132    if (keepalm) {
     133      cout << "--- Writing C_l vector to Output PPF file " << outfile << endl;
     134      ppo.PutObject(clvec, "Cl");
     135      cout << "--- Writing Alm Matrix to Output PPF file " << outfile << endl;
     136      ppo.PutObject(almmtx, "Alm");     
     137    }
     138    else {
     139      cout << "--- Writing C_l vector to Output PPF file " << outfile << endl;
     140      ppo << clvec;
     141    }
    119142  }
    120143}
     
    135158  bool fgfitsout = false;
    136159  bool fgr4 = false;
     160  bool fgkeepalm = false;
    137161  cout << " map2cl : Decoding command line options ... " << endl;
    138162
     
    156180    else if (strcmp(arg[k], "-fitsout") == 0) {
    157181      fgfitsout = true; 
     182    }
     183    else if (strcmp(arg[k], "-keepalm") == 0) {
     184      fgkeepalm = true; 
    158185    }
    159186    else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
     
    174201    if (fgr4) {
    175202      cout << " SphereHEALPix<r_4>  --> Power spectrum C_l<r_4> (float)" << endl;
    176       _Map2Cl<r_4>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder, fgfitsin, fgfitsout);
     203      _Map2Cl<r_4>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder,
     204                              fgfitsin, fgfitsout, fgkeepalm);
    177205    }
    178206    else {
    179207      cout << " SphereHEALPix<r_8>  --> Power spectrum C_l<r_8> (double)" << endl;
    180       _Map2Cl<r_8>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder,fgfitsin, fgfitsout);
     208      _Map2Cl<r_8>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder,
     209                              fgfitsin, fgfitsout, fgkeepalm);
    181210    }
    182211  }
Note: See TracChangeset for help on using the changeset viewer.