Changeset 4027 in Sophya


Ignore:
Timestamp:
Oct 17, 2011, 10:56:41 AM (14 years ago)
Author:
ansari
Message:

Implementatiom prise en compte dA(z) ds la calcul de bruit Pnoise(k) - Reza 17/10/2011

Location:
trunk/Cosmo/RadioBeam
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/pknoise.cc

    r4026 r4027  
    66
    77  Usage:  pknoise pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName
    8           -parname : -noise , -renmax , -scut , -ngen , -z , -bsize , -cszmpc , -prt
     8          -parname : -noise , -renmax, -scut, -ngen, -z, -bsize, -cszmpc, -nbin -prt
    99---------------------------------------------------------------  */
    1010
    … …  
    3030#include "timing.h"
    3131#include "ctimer.h"
     32#include "slininterp.h"
    3233
    3334typedef DR48RandGen RandomGenerator ;
    … …  
    4748       << "       NGen=-1 -> Call ComputeNoisePk(), dont change Four3DPk cell size \n"
    4849       << "    -bsize sx,sy,sz : 3D real space box size (default=512x512x256) \n"   
    49        << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x3) \n"   
    50        << "    -z z,dA,H_z : redshift,ComovDist(z),H(z) (default=1.0,3330,120.5) \n"
     50       << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"   
     51       << "    -z redshift : redshift  (default=1.0) \n"
    5152       << "    -scut SCutValue (default= -100.) \n"
    5253       << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
    53        << "    -prt PrtLev,PrtModulo (default=0,10) " << endl;
     54       << "    -prt PrtLev,PrtModulo (default=0,10) \n"
     55       << "    -nbin PkNBin : nomber of histogram bins for Pnoise(k) (default= 256) \n"
     56       << "  NOTE : this program needs redshift_da.txt andredshift_Hz.txt files" << endl;
    5457  return;
    5558}
    … …  
    8790  double H_z=120.5;  // Hubble_param(z)
    8891  int box3dsz[3]={512,512,256};
    89   double cellsz[3]={3.5,3.5,3.5};
    90 
     92  double cellsz[3]={3.,3.,1.5};
     93  int NBINPK=256;
    9194  int prtlev=0;
    9295  int prtmod=10;
    … …  
    109112    }
    110113    else if (strcmp(arg[ka],"-z")==0) {
    111       sscanf(arg[ka+1],"%lg,%lg,%lg",&z_Redshift,&comov_dA_z,&H_z);  ka+=2;
     114      z_Redshift=atof(arg[ka+1]);  ka+=2;
    112115    }
    113116    else if (strcmp(arg[ka],"-bsize")==0) {
    … …  
    116119    else if (strcmp(arg[ka],"-cszmpc")==0) {
    117120      sscanf(arg[ka+1],"%lg,%lg,%lg",cellsz,cellsz+1,cellsz+2);   ka+=2;
     121    }
     122    else if (strcmp(arg[ka],"-nbin")==0) {
     123      NBINPK=atoi(arg[ka+1]);    ka+=2;
    118124    }
    119125    else if (strcmp(arg[ka],"-prt")==0) {
    … …  
    147153    cout << " pknoise[0] : Executing, output file= " << outfile << endl; 
    148154    POutPersist po(outfile);
     155   
     156    cout << " pknoise[0.b] : Creating interpolating objects from redshift_da.txt and redshift_Hz.txt " << endl;
     157    SLinInterp1D redshift2da;
     158    string interpfilename="redshift_da.txt";
     159    redshift2da.ReadXYFromFile(interpfilename);
     160    SLinInterp1D redshift2hz;
     161    interpfilename="redshift_hz.txt";
     162    redshift2hz.ReadXYFromFile(interpfilename);
     163    comov_dA_z=redshift2da(z_Redshift);
     164    H_z=redshift2hz(z_Redshift);
     165    cout << " Redshift= " << z_Redshift << " -> ComovdA= " << comov_dA_z
     166         << " Mpc  H(z)= " << H_z << " km/s/Mpc " << endl;
    149167
    150168    H21Conversions conv;
    … …  
    184202    RandomGenerator rg;
    185203
     204    double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
     205    double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
     206
    186207    double dkxmpc=2.*DeuxPI;
    187208    double dkympc=2.*DeuxPI;
    188209    double dkzmpc=2.*DeuxPI;
    189210    double angscale=1.;
     211    Vector angscales(box3dsz[2]);
     212    angscales=angscale;
    190213    if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[3]=128; }
    191214    else {
    192215      angscale=comov_dA_z;
     216      for(int kz=0; kz<box3dsz[2]; kz++) angscales(kz)=redshift2da(1420.4/(freq0+(double)kz*dfreq)-1.);
    193217      dkxmpc = DeuxPI/(double)box3dsz[0]/cellsz[0];
    194218      dkympc = DeuxPI/(double)box3dsz[1]/cellsz[1];
    195219      dkzmpc = DeuxPI/(double)box3dsz[2]/cellsz[2];
    196220    }
    197     cout << " pknoise[2.b]: ComovDist" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl;
    198     cout << " pknoise[2.c]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << ","
     221    cout << " pknoise[2.b]: ComovDist=" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl;
     222    cout << " pknoise[2.c]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")="
     223         << conv.getFrequency() << " zmin=" << 1420.4/freq0-1. << " zmax=" << 1420.4/(freq0+box3dsz[2]*dfreq)-1. << endl;
     224    cout << " pknoise[2.d]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << ","
    199225         << box3dsz[2] << ")" << endl;
    200226    Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
    … …  
    203229    m3d.SetPrtLevel(prtlev,prtmod);
    204230
    205     cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    206231    if (NMAX>0) {
     232      cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    207233      PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
    208       double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
    209       double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
    210       cout << " pknoise[3.b]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")="
    211            << conv.getFrequency() << " AngScale=" << angscale << endl;
    212234      pkn.SetFreqRange(freq0, dfreq);
    213       pkn.SetAngScaleConversion(angscale);
     235      pkn.SetAngScaleConversion(angscales);
    214236      pkn.SetPrtLevel(prtlev,prtmod);
    215       HProf hpn = pkn.Compute();
    216       cout << " pknoise[3.c]: writing hpn noise profile to " << outfile << endl;
    217       po << hpn;
     237      HProf hpn = pkn.Compute(NBINPK);
    218238    }
    219239    else {
    220       Histo fracmodok;
    221       DataTable dtnoise;
    222       HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,angscale,SCut);
    223       HProf h1dnoise=arep_p->GetProjNoiseLevel();
    224       HProf h1drep=arep_p->GetProjResponse();
    225       cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
    226       po << PPFNameTag("dtnoise") << dtnoise;
    227       po << PPFNameTag("hpnoise") << hpn;
    228       po << PPFNameTag("fracmodok") << fracmodok;
    229       po << PPFNameTag("h1dnoise") << h1dnoise;
    230       po << PPFNameTag("h1drep") << h1drep;
    231       po << PPFNameTag("h2drep") << h2drep;
    232     }
     240      cout << " pknoise[3]: Computing Noise P(k) m3d.ComputeNoisePk(...) " << endl;
     241      HProf hpn = m3d.ComputeNoisePk(*(arep_p),angscale,SCut,NBINPK);
     242    }
     243       
     244    HProf hpnoise=m3d.GetPk();
     245    DataTable dtnoise;
     246    Histo fracmodok=m3d.FillPkDataTable(dtnoise);
     247    HProf h1dnoise=arep_p->GetProjNoiseLevel();
     248    HProf h1drep=arep_p->GetProjResponse();
     249    cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
     250    po << PPFNameTag("dtnoise") << dtnoise;
     251    po << PPFNameTag("hpnoise") << hpnoise;
     252    po << PPFNameTag("fracmodok") << fracmodok;
     253    po << PPFNameTag("h1dnoise") << h1dnoise;
     254    po << PPFNameTag("h1drep") << h1drep;
     255    po << PPFNameTag("h2drep") << h2drep;
     256 
    233257    rc = 0;
    234258  }  // End of try bloc

  • trunk/Cosmo/RadioBeam/radutil.h

    r3947 r4027  
    5353
    5454  static double SpeedOfLight_Cst;       // Speed of light  m/sec
    55   static double Freq021cm_Cst;       // Speed of light  m/sec
     55  static double Freq021cm_Cst;       // 21 cm transition frequency in MHz 
    5656  static double k_Boltzman_Cst;         // Boltzmann constant  (SI Units)
    5757

  • trunk/Cosmo/RadioBeam/specpk.cc

    r4026 r4027  
    127127  SetPrtLevel();
    128128  SetCellSize();
     129  hp_pk_p_=NULL;  hmcnt_p_=NULL;  hmcntok_p_=NULL;   s2cut_=0.;
    129130}
    130131// Constructor
    … …  
    134135  SetPrtLevel();
    135136  SetCellSize();
    136 }
    137 
     137  hp_pk_p_=NULL;  hmcnt_p_=NULL;  hmcntok_p_=NULL;   s2cut_=0.;
     138}
     139
     140// Destructor
     141Four3DPk::~Four3DPk()
     142{
     143  if (hp_pk_p_) delete hp_pk_p_;
     144  if (hmcnt_p_) delete hmcnt_p_;
     145  if (hmcntok_p_) delete hmcntok_p_;
     146}
    138147
    139148// Generate mass field Fourier Coefficient
    … …  
    200209
    201210// Generate mass field Fourier Coefficient
    202 void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale)
     211void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales)
    203212// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
    204213// typically = ComovRadialDistance
    205214{
     215  if (angscales.Size() != fourAmp.SizeZ())
     216    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
    206217  H21Conversions conv;
    207218  // fourAmp represent 3-D fourier transform of a real input array.
    … …  
    212223    conv.setFrequency(f0+kz*df);
    213224    resp.setLambda(conv.getLambda());
     225    double angsc=angscales(kz);
     226    if (prtlev_>2) 
     227      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc << endl;
    214228    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
    215229      kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
    216230      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
    217231        kxx=(double)kx*dkx_;
    218         rep = resp(kxx*angscale, kyy*angscale);
    219         if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
     232        rep = resp(kxx*angsc, kyy*angsc);
     233        if (rep<1.e-19)  fourAmp(kx, ky, kz) = complex<TF>(9.e19,0.);
    220234        else {
    221235          amp = 1./sqrt(rep)/sqrt(2.);
    … …  
    264278// cells with amp^2=re^2+im^2>s2cut are ignored
    265279// Output : power spectrum (profile histogram)
    266 HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax)
     280HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax, bool fgmodcnt)
    267281{
    268282  // The second half of the array along Y (matrix rows) contain
    … …  
    279293    kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
    280294  }
    281   if (nbin<2) nbin=128;
    282   HProf hp(kmin, kmax, nbin);
    283   hp.SetErrOpt(false);
    284   ComputePkCumul(hp, s2cut);
    285   return hp;
     295  if (nbin<2) nbin=256;
     296  hp_pk_p_ = new HProf(kmin, kmax, nbin);
     297  hp_pk_p_->SetErrOpt(false);
     298  if (fgmodcnt) {
     299    hmcnt_p_ = new Histo(kmin, kmax, nbin);
     300    hmcntok_p_ = new Histo(kmin, kmax, nbin);
     301  }
     302  s2cut_=s2cut;
     303  ComputePkCumul();
     304  return *hp_pk_p_;
    286305}
    287306
    288307// Compute power spectrum as a function of wave number k
    289308// Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
    290 void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
     309void Four3DPk::ComputePkCumul()
    291310{
    292311  uint_8 nmodeok=0;
    … …  
    303322        double kxx=(double)kx*dkx_;
    304323        complex<TF> za = fourAmp(kx, ky, kz);
    305         if (za.real()>8.e9) continue;
     324        //      if (za.real()>8.e19) continue;
    306325        double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
    307326        double amp2 = za.real()*za.real()+za.imag()*za.imag();
    308         if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
    309         hp.Add(wk, amp2);
     327        if (hmcnt_p_) hmcnt_p_->Add(wk);
     328        if ((s2cut_>1.e-9)&&(amp2>s2cut_))  continue;
     329        if (hmcntok_p_) hmcntok_p_->Add(wk);
     330        hp_pk_p_->Add(wk, amp2);
    310331        nmodeok++;
    311332      }
    312333    }
    313334  }
    314   if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
     335  if ((prtlev_>1)||((prtlev_>0)&&(s2cut_>1.e-9))) {
    315336    cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
    316337         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
    … …  
    325346// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
    326347// typically = ComovRadialDistance
    327 HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt,
    328                                double angscale, double s2cut, int nbin, double kmin, double kmax)
     348HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, double angscale, double s2cut, int nbin, double kmin, double kmax)
    329349{
    330350  // The second half of the array along Y (matrix rows) contain
    … …  
    341361    kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
    342362  }
    343   if (nbin<2) nbin=128;
    344   HProf hp(kmin, kmax, nbin);
    345   hp.SetErrOpt(false);
    346   Histo hmcnt(kmin, kmax, nbin);
    347   Histo hmcntok(kmin, kmax, nbin);
     363  if (nbin<2) nbin=256;
     364  hp_pk_p_ = new HProf(kmin, kmax, nbin);
     365  hp_pk_p_->SetErrOpt(false);
     366  hmcnt_p_ = new Histo(kmin, kmax, nbin);
     367  hmcntok_p_ = new Histo(kmin, kmax, nbin);
     368  s2cut_=s2cut;
    348369
    349370  uint_8 nmodeok=0;
    … …  
    362383        double rep=resp(kxx*angscale, kyy*angscale);
    363384        double amp2 = (rep>1.e-19)?1./rep:1.e19;
    364         hmcnt.Add(wk);
    365         if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
    366         hmcntok.Add(wk);
    367         hp.Add(wk, amp2);
     385        hmcnt_p_->Add(wk);
     386        if ((s2cut_>1.e-9)&&(amp2>s2cut_))  continue;
     387        hmcntok_p_->Add(wk);
     388        hp_pk_p_->Add(wk, amp2);
    368389        nmodeok++;
    369390      }
    … …  
    375396  }
    376397
    377   fracmodok=hmcntok/hmcnt;
    378 
     398  return *hp_pk_p_;
     399}
     400
     401// Fills a data table from the computed P(k) profile histogram and mode count
     402Histo Four3DPk::FillPkDataTable(DataTable& dt)
     403{
     404  if (hp_pk_p_==NULL) throw ParmError("Four3DPk::FillPkDataTable P(k) NOT computed");
     405  if ((hmcnt_p_==NULL)||(hmcntok_p_==NULL)) throw ParmError("Four3DPk::FillPkDataTable Mode count histos NOT filled");
     406  HProf& hp=(*hp_pk_p_);
     407  Histo& hmcnt=(*hmcnt_p_);
     408  Histo& hmcntok=(*hmcntok_p_);
     409  Histo fracmodok=hmcntok/hmcnt;
    379410  char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
    380411  dt.Clear();
    … …  
    393424    dt.AddRow(dtr);
    394425  }
    395   return hp;
     426  return fracmodok;
    396427}
    397428
    … …  
    403434PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen,
    404435                                     const char* tit)
    405   : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
     436  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ())
    406437{
    407438  SetFreqRange();
    … …  
    410441}
    411442
    412 HProf PkNoiseCalculator::Compute()
     443HProf PkNoiseCalculator::Compute(int nbin, double kmin, double kmax)
    413444{
    414445  Timer tm(title.c_str());
    … …  
    416447  HProf hnd;
    417448  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT 
    418        << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscale=" << angscale_ << endl;
     449       << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscales=" << angscales_(0)
     450       << " ... " << angscales_(angscales_.Size()-1) << endl;
    419451  for(int igen=0; igen<NGEN; igen++) {
    420     pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscale_);
    421     if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
    422     else pkn3d.ComputePkCumul(hnd,S2CUT);
     452    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_);
     453    if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
     454    else pkn3d.ComputePkCumul();
    423455    if ((prtlev_>0)&&(igen>0)&&(((igen-1)%prtmodulo_)==0))
    424456      cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
    425457  }
    426   return hnd;
     458  return pkn3d.GetPk() ;
    427459}
    428460

  • trunk/Cosmo/RadioBeam/specpk.h

    r4026 r4027  
    4646  Four3DPk(TArray< complex<TF> > & fourcoedd, RandomGeneratorInterface& rg);
    4747  Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx=128, sa_size_t szy=256, sa_size_t szz=128);
     48  virtual ~Four3DPk();
     49
    4850  inline void SetCellSize(double dkx=DeuxPI, double dky=DeuxPI, double dkz=DeuxPI)
    4951  { dkx_=dkx;  dky_=dky;  dkz_=dkz; }
    … …  
    5456// typically = ComovRadialDistance
    5557  void ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale=1., bool crmask=false);
    56   void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale=1.);
     58  void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales);
    5759
    5860// Return the array size
    5961  inline sa_size_t NCells() { return fourAmp.Size(); }
     62  inline sa_size_t SizeX() { return fourAmp.SizeX(); }
     63  inline sa_size_t SizeY() { return fourAmp.SizeY(); }
     64  inline sa_size_t SizeZ() { return fourAmp.SizeZ(); }
     65
    6066// Set the cell size/step in Fourier Space
    6167// Return the fourier amplitude matrix 
    … …  
    6672
    6773// Return the reconstructed power spectrum as a profile histogram   
    68   HProf ComputePk(double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
    69   void  ComputePkCumul(HProf& hp, double s2cut=0.);
     74  HProf ComputePk(double s2cut=0., int nbin=256, double kmin=0., double kmax=-1., bool fgmodcnt=false);
     75  void  ComputePkCumul();
    7076
    7177// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
    7278// typically = ComovRadialDistance
    73   HProf ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt, double angscale=1.,
    74                        double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
     79  HProf ComputeNoisePk(Four2DResponse& resp, double angscale=1., double s2cut=0.,
     80                       int nbin=256, double kmin=0., double kmax=-1.);
     81
     82  // Fills a data table from the computed P(k) profile histogram and mode count
     83  Histo FillPkDataTable(DataTable& dt);
     84  inline HProf& GetPk() { return *hp_pk_p_; }
    7585
    7686protected:
    … …  
    8191  int prtlev_;
    8292  int prtmodulo_;
     93  // Profile histograms for power spectrum and number of modes
     94  HProf* hp_pk_p_;
     95  Histo* hmcnt_p_;
     96  Histo* hmcntok_p_;
     97  double s2cut_;
    8398};
    8499
    … …  
    94109  { freq0_=freq0;  dfreq_=dfreq; }
    95110  inline void SetAngScaleConversion(double angscale=1.)
    96   { angscale_=angscale; }
     111  { angscales_=angscale; }
     112  inline void SetAngScaleConversion(Vector& angscs)
     113  { angscales_=angscs; }
    97114  inline void SetS2Cut(double s2cut=100.)
    98115  {  S2CUT=s2cut; }
    99116  inline double GetS2Cut() { return S2CUT; }
    100   HProf Compute();
     117  HProf Compute(int nbin=256, double kmin=0., double kmax=-1.);
    101118  inline int SetPrtLevel(int lev=0, int prtmod=10)
    102119  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
    … …  
    106123  Four2DResponse& frep;
    107124  double freq0_,dfreq_;
    108   double angscale_;
     125  Vector angscales_;
    109126  double S2CUT;
    110127  int NGEN;
Note: See TracChangeset for help on using the changeset viewer.