// Implementation of all routines concerning a model photon. The
// object 'photon follows a photon's life from its creation in some
// emission region of the model space, through scattering processes,
// until either its absorption or detection. All necessary routines
// and parameters are defined here, including the procedures for
// electron scattering or scattering off a generalized dust
// grain. Some features are marked with a symbol (*) denoting that
// they are not used by the program at its current state but that they
// are reserved for future purposes or diagnostic reasons.
//
// STOKES, version 1.0, Nov 2004


#include <stdlib.h>
#include <math.h>
#include "photon-stokes-v1.0.h"


// set the photon parameters position, flying direction, Stokes vector, weight,
// and wavelength to the given values (*)
void Photon::SetPhotonData(double setpos[],double setdirect[][3],
                           double setI,double setQ,double setU,double setV,
                           double setW,double setlambda)
{
  int i;
  Equal(pos,setpos);
  for(i = 0; i <= 2; i++)
  {
    Equal(direct[i],setdirect[i]);
  }
  I = setI;
  Q = setQ;
  U = setU;
  V = setV;
  W = setW;
  lambda = setlambda;
  abs = false;
}


// reset photon to an initial state, located in a continuum or line
// emission region, being unpolarized, and flying toward a random
// direction
void Photon::Reset(Model M)
{

  // origin of the model space
  const double origin[3] = {0,0,0};

  // default frame of the model space, xyz coordinate system
  const double frame[3][3] = {{1,0,0},{0,1,0},{0,0,1}};

  int i,j,line;
  double R,phi,theta,alpha,min,max,vr[3],r;
  CTRegion CTR;
  BLRegion BLR;
  NLRegion NLR;

  // set the co-moving frame of the photon equal to the
  // default coordinate system of the model space
  for(i=0; i<=2; i++)
  {
    for(j=0; j<=2; j++)
    {
      direct[i][j] = frame[i][j];
    }
  }

  // sample phi and theta for the flight direction assuming isotropic emission
  phi = double(rand())/RAND_MAX * 2*pi;
  do
  {
    theta = acos(1 - 2*double(rand())/RAND_MAX);
  }
  while ((theta == 0) || (theta == pi));

  // turn the photon's frame into the flight direction
  TurnDirection(2,phi);
  TurnDirection(1,theta-pi/2);

  // decide if the photon to create is a continuum photon, a broad
  // line photon, or a narrow line photon
  CTR = M.GetCTR();
  BLR = M.GetBLR();
  NLR = M.GetNLR();
  r = double(rand())/RAND_MAX;
  if (r < NLR.rli)
  {
    line = 2;
  }
  else
  {
    if (r < NLR.rli+BLR.rli)
    {
      line = 1;
    }
    else
    {
      line = 0;
    }
  }
  
  switch (line)
  {

    // continuum photon, sample starting position inside the
    // cylinder CTR
    case 0:  
      if (CTR.R==0)
      {
        for(i=0; i<=2; i++)
        {
          pos[i] = origin[i];
        }
      }
      else
      {
        r = double(rand())/RAND_MAX;
        R = sqrt((CTR.R-CTR.b)*(CTR.R-CTR.b)+4*r*CTR.R*CTR.b);
        r = double(rand())/RAND_MAX;
        phi = 2*pi*r;
        pos[0] = R*cos(phi);
        pos[1] = R*sin(phi);
        r = double(rand())/RAND_MAX;
        pos[2] = CTR.a*(2*r-1);
      }

      // sample the photon wavelength according to an I ~ nu^(-alpha)
      // intensity spectrum; shift the wavelength if the emission
      // region is in motion
      r = double(rand())/RAND_MAX;
      if (CTR.alpha!=0)
      {
        min = pow(M.LambdaMin(),CTR.alpha);
        max = pow(M.LambdaMax(),CTR.alpha);
        lambda = pow(min+r*(max-min),1/CTR.alpha);
      }
      else
      {
        min = M.LambdaMin();
        max = M.LambdaMax();
        lambda = min*pow(max/min,r);
      }
      if (CTR.vr!=0)
      {
        SMult(vr,CTR.vr/GetDistance(),pos);
        LambdaShift(false,vr);
      }
      break;

    // broad line photon, sample starting position inside the
    // cylinder BLR and wavelength according to a lorentz profile;
    // shift the wavelength if the emission region is in motion
    case 1:  
      r = double(rand())/RAND_MAX;
      R = sqrt((BLR.R-BLR.b)*(BLR.R-BLR.b)+4*r*BLR.R*BLR.b);
      r = double(rand())/RAND_MAX;
      phi = 2*pi*r;
      pos[0] = R*cos(phi);
      pos[1] = R*sin(phi);
      r = double(rand())/RAND_MAX;
      pos[2] = BLR.a*(2*r-1);
      lambda = GetLineLambda(M,1,r);
      if (BLR.vr!=0)
      {
        SMult(vr,BLR.vr/GetDistance(),pos);
        LambdaShift(false,vr);
      }
      break;

    // narrow line photon, sample starting position inside the
    // double-cone NLR and wavelength according to a lorentz profile;
    // shift the wavelength if the emission region is in motion
    default:
      r = double(rand())/RAND_MAX;
      R = pow(r*pow(NLR.b,3)+(1-r)*pow(NLR.a,3),1.0/3);
      r = double(rand())/RAND_MAX;
      phi = 2*pi*r;
      r = double(rand())/RAND_MAX;
      theta = acos(1-r*(1-cos(NLR.theta)));
      CartCoord(pos,phi,theta);
      r = double(rand())/RAND_MAX;
      if (r<0.5)
      {
        SMult(pos,R,pos);
      }
      else
      {
        SMult(pos,-R,pos);
      }
      lambda = GetLineLambda(M,2,r);
      if (NLR.vr!=0)
      {
        SMult(vr,NLR.vr/GetDistance(),pos);
        LambdaShift(false,vr);
      }
      break;
  }

  // initialize the Stokes vector to an unpolarized photon
  I = 1;
  V = 0;
  Q = 0;
  U = 0;

  // the weight factor may be used in the future, currently it
  // is set to unity for all photons
  W = 1;

  // the photon is outside any scattering region
  ScattFlag = ScattMax;

  // no scatterings occurred yet
  ScattNum = 0;

  // assume electron scattering cross-section
  CS = TCS;

  // the photon is not absorbed
  abs = false;
}


// sample a wavelength for a line photon, decide according to the
// parameter 'line' if the photon is a broad or a narrow line
// photon. The wavelength is sampled assuming a Lorentz-profile of the
// intensity versus frequency spectrum
double Photon::GetLineLambda(Model M,int line,double r)
{
  double lambda,l1,l2,l,f;
  BLRegion BLR;
  NLRegion NLR;

  // broad line photon
  if (line == 1)
  {
    BLR = M.GetBLR();
    l1 = M.LambdaMin();
    l2 = M.LambdaMax();
    do
    {
      l = (l1+l2)/2;
      f = 0.5*log(pow((l-BLR.lambda),2)+pow(BLR.width,2))+
          BLR.lambda/BLR.width*atang(l-BLR.lambda,BLR.width)-
          BLR.low-r*BLR.norm;
      if (f >= 0)
        l2 = l;
      else
        l1 = l;
    }
    while (l2-l1>LinePrecision);
    lambda = (l1+l2)/2;
  }

  // narrow line photon
  else
  {
    NLR = M.GetNLR();
    l1 = M.LambdaMin();
    l2 = M.LambdaMax();
    do
    {
      l = (l1+l2)/2;
      f = 0.5*log(pow((l-NLR.lambda),2)+pow(NLR.width,2))+
          NLR.lambda/NLR.width*atang(l-NLR.lambda,NLR.width)-
          NLR.low-r*NLR.norm;
      if (f >= 0)
        l2 = l;
      else
        l1 = l;
    }
    while (l2-l1>LinePrecision);
    lambda = (l1+l2)/2;
  }

  return lambda;
}


// if the scattering region Cloud is in the photon's flight direction,
// then set 'tpre' to the distance to the cloud's closest border,
// otherwise leave 'tpre' unchanged
void Photon::FindNextBorder(scatterer Cloud,double& tpre,double pos[],
                            double direct[])
{
  double t;

  t = tpre;
  switch(Cloud.shape)
  {

    // check for a toroidal scattering region
    case 0:
      TorusCoinc(pos,direct,t,Cloud.a,Cloud.b,Cloud.R);
      break;

    // check for a cylindrical scattering region, consisting of four
    // connected surfaces
    case 1:
      CylinderCoinc(pos,direct,t,Cloud.R-Cloud.b,Cloud.a);
      CylinderCoinc(pos,direct,t,Cloud.R+Cloud.b,Cloud.a);
      RingCoinc(pos,direct,t,Cloud.R-Cloud.b,Cloud.R+Cloud.b,-Cloud.a);
      RingCoinc(pos,direct,t,Cloud.R-Cloud.b,Cloud.R+Cloud.b,Cloud.a);
      break;

    // check for a double-conical scattering region, consisting of three
    // connected surfaces
    case 2:
      ConeCoinc(pos,direct,Cloud.a,Cloud.b,Cloud.theta,t);
      SphereSegCoinc(pos,direct,Cloud.pos,t,Cloud.a,Cloud.theta);
      SphereSegCoinc(pos,direct,Cloud.pos,t,Cloud.b,Cloud.theta);
      break;

    // check for a spherical scattering region
    default:
      SphereSegCoinc(pos,direct,Cloud.pos,t,Cloud.a,pi/2);
      break;

  }

  // find the closest scattering region
  if (t<tpre)
  {
    tpre = t;
  }

}


// search for the nearest scattering region along the photons flight direction,
// if there is such cloud, then move the photon to the clouds border
bool Photon::ShiftToNextScatt(Model M,Mie Dust)
{
  int i;
  double t,tpre,NewPos[3];
  scatterer Cloud;

  // set maximum flight distance that is possible
  tpre = M.GetMaxExp();
  t = M.GetMaxExp();

  // assume photon to be inside no scattering region
  ScattFlag = ScattMax;

  // loop over all defined scattering regions
  for(i=0; i<M.GetScattNum(); i++)
  {
    Cloud = M.GetCloud(i);

    // if the photon already is inside some scattering cloud
    if (InsideCloud(Cloud))
    {
      ScattFlag = i;

      // do not displace the photon
      t = -M.GetGap();

      break;
    }

    // search for the closest scattering region in the photons
    // flight direction
    else
    {
      FindNextBorder(Cloud,t,pos,direct[0]);

      // if the distance to the Cloud is smaller than the previous one
      // then re-define the distance
      if (t<tpre)
      {
        ScattFlag = i;
        tpre = t;
      }
    }
  }

  // new scattering region found
  if (ScattFlag < ScattMax)
  {

    // chose the scattering cross-section (unit ccm/parsec)
    Cloud = M.GetCloud(ScattFlag);
    if (Cloud.dust)
    {
      CS = 3.08633e18*Dust.GetCext(lambda);
    }
    else
    {
      CS = TCS;
    }

    // displace the photon, add 'Gap' to make sure that the photon is
    // declared to be inside the cloud
    SMult(NewPos,t+M.GetGap(),direct[0]);
    VecSum(NewPos,pos,NewPos);
    SetPosition(NewPos);

    return true;
  }

  // no scattering region found
  else
  {
    return false;
  }

}


// sample a path length for the photon and shift it within the given
// scattering region. If the photon leaves the border of the region
// return false, otherwise return true
bool Photon::Proceed(Model M)
{
  double s,r,t,density,connect[3],OldPos[3];
  scatterer Cloud;

  // save the old photon position
  Equal(OldPos,pos);

  // sample a path length within the current scattering region
  Cloud = M.GetCloud(ScattFlag);
  density = Cloud.Density;
  do
  {
    r = double(rand())/RAND_MAX;
  }
  while (r==0);
  s = -log(r)/(density*CS);

  // displace the photon to the new position
  SMult(connect,s,direct[0]);
  VecSum(pos,pos,connect);

  // photon stays inside the scattering region
  if (InsideCloud(Cloud))
  {
    return true;
  }

  // photon leaves the scattering region
  else
  {

    // return to the old photon position
    SetPosition(OldPos);

    // shift just outside the border of the current
    // scattering region
    t = s;
    FindNextBorder(Cloud,t,pos,direct[0]);
    SMult(connect,t+M.GetGap(),direct[0]);
    VecSum(pos,pos,connect);

    return false;
  }
}


// return true if the photon is located inside the scattering region Cloud,
// return false otherwise;
bool Photon::InsideCloud(scatterer Cloud)
{
  double a,b,con[3],r,tn;
  bool inside;

  // assume photon outside cloud
  inside = false;
  switch(Cloud.shape)
  {

    // check for a toroidal scattering region
    case 0:
      a = Cloud.a;
      b = Cloud.b;
      if (pos[2]*pos[2]/(a*a)+pow(Cloud.R-sqrt(pos[0]*pos[0]+pos[1]*pos[1]),2)/
         (b*b)<=1)
      {
        inside = true;
      }
      break;

    // check for a cylindrical scattering region
    case 1:
      r = sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
      if ((r>=Cloud.R-Cloud.b)&&(r<=Cloud.R+Cloud.b)&&(fabs(pos[2])<=Cloud.a))
      {
        inside = true;
      }
      break;

    // check for a double-conical scattering region
    case 2:
      r = SProd(pos,pos);
      tn = tan(pi/2-Cloud.theta);
      if ((tn*tn*(pos[0]*pos[0]+pos[1]*pos[1])<pos[2]*pos[2])&&
          (r>=Cloud.a*Cloud.a)&&(r<=Cloud.b*Cloud.b))
      {
        inside = true;
      }
      break;

    // check for a spherical scattering region
    default:
      VecDiff(con,pos,Cloud.pos);
      if (SProd(con,con) <= Cloud.a*Cloud.a)
      {
        inside = true;
      }
      break;

  }
  return inside;
}


// turn the photon's flight direction by the angle 'alpha' and around the
// coordinate axis given by int; mathematically this is a matrix multiplication
// of the rotation matrix mat with each unity vector direct[i] of the photon's
// co-moving coordinate system
void Photon::TurnDirection(int axis,double ep)
{
  int i,j;
  double sn,cs,temp,norm;
  double tempdirect[3];
  double mat[3][3];

  sn = sin(ep);
  cs = cos(ep);

  // prepare the rotation matrix mat
  for(i = 0; i <=2; i++) mat[i][i] = pow(direct[axis][i],2) * (1 - cs) + cs;
  temp = direct[axis][0] * direct[axis][1] * (1 - cs);
  mat[0][1] = temp - direct[axis][2] * sn;
  mat[1][0] = temp + direct[axis][2] * sn;
  temp = direct[axis][0] * direct[axis][2] * (1 - cs);
  mat[0][2] = temp + direct[axis][1] * sn;
  mat[2][0] = temp - direct[axis][1] * sn;
  temp = direct[axis][1] * direct[axis][2] * (1 - cs);
  mat[1][2] = temp - direct[axis][0] * sn;
  mat[2][1] = temp + direct[axis][0] * sn;

  // matrix multiplication and re-normalization of the unity vectors
  // to avoid numerical inaccuracies
  for(i = 0; i <= 2; i++)
  {
    MatVecMult(mat,direct[i],tempdirect);
    for(j = 0; j <= 2; j++) direct[i][j] = tempdirect[j];
    norm = sqrt(SProd(direct[i],direct[i]));
    for(j = 0; j <= 2; j++) direct[i][j] = 1/norm * direct[i][j];
  }

}


// turn the photon's polar plane by the angle 'alpha' around the axis given by
// its flight direction; this is a rotation of the polarization plane, defined
// by direct[1] and direct[2], within itself; the Stokes vector is adjusted to
// keep the correct polarization
void Photon::TurnPolarPlane(double alpha)
{
  double QNew,UNew,c1,c2;

  c1 = cos(2*alpha);
  c2 = sin(2*alpha);
  QNew = c1*Q - c2*U;
  UNew = c2*Q + c1*U;
  Q = QNew;
  U = UNew;
  TurnDirection(0,alpha);
}


// perform an electron scattering event and change all affected
// photon parameters
void Photon::ElectronScattering(Model M)
{
  double PL,theta,phi,vr[3];
  scatterer Cloud;

  Cloud = M.GetCloud(ScattFlag);

  // if the scattering cloud is in motion transform
  // the photon wavelength
  if (Cloud.vr!=0)
  {
    SMult(vr,Cloud.vr/GetDistance(),pos);
    LambdaShift(true,vr);
  }

  // sample a scattering angle theta
  theta = GetElecScattTheta(double(rand())/RAND_MAX);

  // compute the incident polarization and sample a scattering angle phi
  PL = sqrt(Q*Q+U*U)/I;
  phi = GetElecScattPhi(double(rand())/RAND_MAX,PL,theta);

  // adjust the photon's flight direction and Stokes parameters
  TurnPolarPlane(pi-0.5*atang(U,Q)+phi);
  ElectTurnTheta(theta);

  // if the scattering cloud is in motion transform
  // the photon wavelength
  if (Cloud.vr!=0)
  {
    LambdaShift(false,vr);
  }

  // increase the number of scatterings
  ScattNum = ScattNum + 1;

}


// sample an electron scattering angle theta according to the
// random number r; the necessary trigonometric equation is
// solved by an iterative method
double  Photon::GetElecScattTheta(double r)
{
  double p1,p2,p,f;

  p1 = 0;
  p2 = pi;
  do
  {
    p = (p2+p1)/2;
    f = 0.5+0.125*cos(p)*(pow(sin(p),2)-4) - r;
    if (f >= 0)
      p2 = p;
    else
      p1 = p;
  }
  while (p2-p1 > pi/180*PFPrecision);
  return (p2+p1)/2;
}


// sample an electron scattering angle phi according to the
// random number r, incident polarization PL, and the scattering
// angle theta; the necessary trigonometric equation is solved by
// an iterative method
double Photon::GetElecScattPhi(double r,double PL,double theta)
{
  double p1,p2,p,f,N,csq;

  p1 = 0;
  p2 = 2*pi;
  N = 1/(2*pi);
  csq = pow(cos(theta),2);
  do
  {
    p = (p2+p1)/2;
    f = N*(p - (1-csq)/(1+csq)*PL*sin(2*p)/2) - r;
    if (f >= 0)
      p2 = p;
    else
      p1 = p;
  }
  while (p2-p1 > pi/180*PFPrecision);
  return (p2+p1)/2;
}


// turn the flight direction of the photon due to electron scattering by
// the angle theta, adjust the Stokes parameters
void Photon::ElectTurnTheta(double theta)
{
  double INew,VNew,QNew,UNew,cs;

  cs = cos(theta);

  // compute new Stokes parameters an normalize them to I=1
  INew = 0.5*(1+cs*cs)*I + 0.5*(cs*cs-1)*Q;
  QNew = 0.5*(cs*cs-1)*I + 0.5*(1+cs*cs)*Q;
  UNew = cs * U;
  VNew = cs * V;
  if (INew!=0)
  {
    Q = QNew/INew;
    U = UNew/INew;
    V = VNew/INew;
    I = 1;
  }

  // turn the flight direction by the scattering angle theta
  TurnDirection(2,theta);

}


// perform a dust scattering event and change all affected
// photon parameters
void Photon::DustScattering(Model M,Mie Dust)
{
  int i;
  double PL,theta,phi,vr[3],r;
  scatterer Cloud;

  Cloud = M.GetCloud(ScattFlag);

  // decide if the photon is absorbed or scattered
  r = double(rand())/RAND_MAX;

  // photon is scattered
  if(r < Dust.GetAlbedo(lambda))
  {

    // if the scattering cloud is in motion transform
    // the photon wavelength
    if (Cloud.vr!=0)
    {
      SMult(vr,Cloud.vr/GetDistance(),pos);
      LambdaShift(true,vr);
    }

    // sample a scattering angle theta
    r = double(rand())/RAND_MAX;
    theta = GetDustScattTheta(r,Dust);

    // compute incident polarization and sample a scattering angle phi
    PL = sqrt(Q*Q+U*U)/I;
    r = double(rand())/RAND_MAX;
    phi = GetDustScattPhi(r,Dust,PL,theta);

    // adjust the photon's flight direction and Stokes parameter
    TurnPolarPlane(pi-0.5*atang(U,Q)+phi);
    DustTurnTheta(Dust,theta);

    // if the scattering cloud is in motion transform
    // the photon wavelength
    if (Cloud.vr!=0)
    {
      LambdaShift(false,vr);
    }

    // increase the number of scatterings
    ScattNum = ScattNum + 1;
  }

  // photon is absorbed
  else
  {
    abs = true;
  }
}


// find a dust scattering angle theta according to the
// random number r, the theta angles are stored in a
// pre-computed look-up table
double Photon::GetDustScattTheta(double r,Mie Dust)
{
  return Dust.GetTheta(lambda,r);
}


// sample a dust scattering angle phi according to the
// random number r, incident polarization PL, and the scattering
// angle theta; the necessary trigonometric equation is solved by
// an iterative method
double Photon::GetDustScattPhi(double r,Mie Dust,double PL,double theta)
{
  double p1,p2,p,f,N;

  p1 = 0;
  p2 = 2*pi;
  N = 1/(2*pi);
  do
  {
    p = (p2+p1)/2;
    f = N*(p-Dust.GetS12(lambda,theta)/Dust.GetS11(lambda,theta)*
           PL*sin(2*p)/2)-r;
    if (f >= 0)
      p2 = p;
    else
      p1 = p;
  }
  while (p2-p1 > pi/180*PFPrecision);
  return (p2+p1)/2;
}


// turn the flight direction of the photon due to dust scattering by
// the angle theta, adjust the Stokes parameters
void Photon::DustTurnTheta(Mie Dust,double theta)
{
  double S11,S12,S33,S34,INew,VNew,QNew,UNew;

  // get the pre-computed Mueller matrix for this case
  Dust.GetMueller(lambda,theta,S11,S12,S33,S34);

  // compute new Stokes parameters an normalize them to I=1
  INew = S11*I+S12*Q;
  QNew = S12*I+S11*Q;
  UNew = S33*U+S34*V;
  VNew = -S34*U+S33*V;
  if (INew!=0)
  {
    Q = QNew/INew;
    U = UNew/INew;
    V = VNew/INew;
    I = 1;
  }

  // turn the flight direction by the scattering angle theta
  TurnDirection(2,theta);
}


// Transform the photon wavelength to (ToFrame=true) or from (ToFrame=false)
// a reference frame moving in the radial direction with the velocity vr;
// the new wavelength is computed by a relativistic Doppler shift formula.
void Photon::LambdaShift(bool ToFrame,double vr[])
{
  double betak,beta2;

  // project on the radial component of the velocity
  // and  compute beta^2
  betak = SProd(vr,direct[0])/vc;
  beta2 = SProd(vr,vr)/(vc*vc);

  // compute the wavelength shift
  if (ToFrame)
  {
    lambda = lambda*(1+betak)/sqrt(1-beta2);
  }
  else
  {
    lambda = lambda*(1-betak)/sqrt(1-beta2);
  }
}


// adjust the Stokes parameters of the photon to the polarization plane of
// a detector in the direction given by theta and phi
void Photon::AdjustFrame(double phi,double theta)
{
  double vp[3],sp,angle,FDirect[3];

  // find the unity vector FDirect of the detector's polarization plane
  // that points 'down-wards'
  if (theta <= pi/2)
  {
    CartCoord(FDirect,phi,theta+pi/2);
  }
  else
  {
    CartCoord(FDirect,phi+pi,3*pi/2-theta);
  }

  // compute the angle between FDirect and the co-moving unity vector
  // direct[2]; from their scalar product sp; check if sp is within the
  // allowed limits
  sp = SProd(FDirect,direct[2]);
  if (sp > 1)
  {
    sp = 1;
  }
  if (sp < -1)
  {
    sp = -1;
  }
  angle = acos(sp);

  // compute the vector product between FDirect and direct[2] to
  // decide in which direction the photon's polarization plane has
  // to be turned
  vp[0] = FDirect[1]*direct[2][2]-direct[2][1]*FDirect[2];
  vp[1] = FDirect[2]*direct[2][0]-direct[2][2]*FDirect[0];
  vp[2] = FDirect[0]*direct[2][1]-direct[2][0]*FDirect[1];
  if (SProd(vp,direct[0]) >= 0)
  {
    TurnPolarPlane(-angle);
  }
  else
  {
    TurnPolarPlane(angle);
  }
}


// return the photon wavelength
double Photon::GetLambda()
{
  return lambda;
}


// set photon position to 'setpos'
void Photon::SetPosition(double setpos[])
{
  int i;
  for(i = 0; i <= 2; i++)
    pos[i] = setpos[i];
}


// set 'getpos' to the current photon position
void Photon::GetPosition(double getpos[])
{
  int i;
  for(i = 0; i <= 2; i++)
    getpos[i] = pos[i];
}


// return the x-component of the flight direction
double Photon::GetDirect1()
{
  return direct[0][0];
}


// return the y-component of the flight direction
double Photon::GetDirect2()
{
  return direct[0][1];
}


// return the z-component of the flight direction
double Photon::GetDirect3()
{
  return direct[0][2];
}


// set the scattering region number to 'NewFlag'
void Photon::SetScattFlag(int NewFlag)
{
  ScattFlag = NewFlag;
}


// return the current scattering region number
int Photon::GetFlag()
{
  return ScattFlag;
}


// return the number of scatterings the photon has undergone (*)
int Photon::GetScattNum()
{
  return ScattNum;
}


// return the distance of the photon from the origin of the model space
double Photon::GetDistance()
{
  return sqrt(SProd(pos,pos));
}


// return the I value of the Stokes vector
double Photon::GetI()
{
  return I;
}


// return the Q value of the Stokes vector
double Photon::GetQ()
{
  return Q;
}


// return the U value of the Stokes vector
double Photon::GetU()
{
  return U;
}


// return the V value of the Stokes vector
double Photon::GetV()
{
  return V;
}


// return the weight factor of the photon (*)
double Photon::GetW()
{
  return W;
}


// return the value of the photon's absorption flag
bool Photon::Absorbed()
{
  return abs;
}