source: Sophya/trunk/Cosmo/RadioBeam/brsrc.cc@ 3803

#include "brsrc.h"
#include "srandgen.h"
#include "fioarr.h"
#include "vector3d.h"

//=================================================
BRSourceGen::BRSourceGen()
{
  int ns = 3*9*5;
  freq.SetSize(ns);
  amp.SetSize(ns);
  angX.SetSize(ns);
  angY.SetSize(ns);
  phase.SetSize(ns);

  amp = 1.;
  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);

  double fr[3] = {0.1, 0.25, 0.4};   // attention 0 <= f < 0.5
  int k = 0;
  for(int kf=0; kf<3; kf++) {
    for(int ja=-4; ja<=4; ja++) {
      for(int ia=-2; ia<=2; ia++) {
        freq(k) = fr[kf];
        // Pour lambda=50 cm, angle M_PI/12. -> delta phi = 2 pi sur ~ 1.9 m
        angX(k) = M_PI/12.*(double)ja; 
        // Pour lambda=50 cm, angle M_PI/5000. -> delta phi = 2 pi sur ~ 800 m
        angY(k) = M_PI/3000.*(double)ia;
        k++;
      }
    }
  }
}

BRSourceGen::BRSourceGen(int ns, double maxangXrad, double maxangYrad, 
                         double minfreq, double minamp, double maxamp)
{
  freq.SetSize(ns);
  amp.SetSize(ns);
  angX.SetSize(ns);
  angY.SetSize(ns);
  phase.SetSize(ns);

  if (minfreq>=0.5) {
    cout << "BRSourceGen::BRSourceGen()/Warning minfreq=" << minfreq << " >=0.5 --> =0." << endl;
    minfreq = 0.;
  }
  for(int i=0; i<ns; i++) {
    freq(i) = minfreq+drand01()*(0.5-minfreq);
    amp(i) = minamp+drand01()*(maxamp-minamp);
    angX(i) = drandpm1()*maxangXrad;
    angY(i) = drandpm1()*maxangYrad;
  }
  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);
}

BRSourceGen::BRSourceGen(vector<double> f, int nsf, double maxangXrad, double maxangYrad, 
                         double minfreq, double minamp, double maxamp)
{
  int ns = nsf*f.size();

  freq.SetSize(ns);
  amp.SetSize(ns);
  angX.SetSize(ns);
  angY.SetSize(ns);
  phase.SetSize(ns);
  if (minfreq>=0.5) {
    cout << "BRSourceGen::BRSourceGen()/Warning minfreq=" << minfreq << " >=0.5 --> =0." << endl;
    minfreq = 0.;
  }

  int j = 0;
  for(int k=0; k<(int)f.size(); k++) {
    if (f[k] >= 0.5) {
      cout << "BRSourceGen::BRSourceGen() f[k=" << k << "]=" << f[k] << " >=0.5 ignored " << endl;
      continue;
    }
    for(int i=0; i<nsf; i++) {
      freq(j) = f[k];
      amp(j) = minamp+drand01()*(maxamp-minamp);
      angX(j) = drandpm1()*maxangXrad;
      angY(j) = drandpm1()*maxangYrad;
      j++;
    }
  }
  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);
}

BRSourceGen::BRSourceGen(string ppfname)
{
  ReadPPF(ppfname);
}

void BRSourceGen::WritePPF(string ppfname)
{
  POutPersist po(ppfname);
  po << freq << amp << angX << angY << phase;
}

void BRSourceGen::ReadPPF(string ppfname)
{
  PInPersist pi(ppfname);
  pi >> freq >> amp >> angX >> angY >> phase;
}

void BRSourceGen::Print(ostream& os)
{
  os << "---- BRSourceGen::Print() - NSources= " << freq.Size() << endl; 
  for(int k=0; k<freq.Size(); k++) {
    os << k << " - f=" << freq(k) << " amp=" << amp(k) << " angX=" << angX(k) 
       << " rad =" << Angle(angX(k)).ToDegree() << " deg angY="
       << " rad =" << Angle(angX(k)).ToArcMin() << " arcmin" 
       << " phase=" << phase(k) << endl; 
  }
}


NTuple BRSourceGen::Convert2Table(double freq0)
{
  const char* nomsntsrc[9] = {"num","freq","freqred","amp","angrad","angdeg","phase",
                              "angyrad","angymin"};
  NTuple nt(9, nomsntsrc);
  double xnt[15];
  for(int kk=0; kk<freq.Size(); kk++) {
    xnt[0] = kk;
    xnt[1] = freq(kk)+freq0;
    xnt[2] = freq(kk);
    xnt[3] = amp(kk);
    xnt[4] = angX(kk);
    xnt[5] = Angle(angX(kk)).ToDegree();
    xnt[6] = phase(kk);
    xnt[7] = angY(kk);
    xnt[8] = Angle(angY(kk)).ToArcMin();
    nt.Fill(xnt);
  }

  return nt;
}