// Implementation of the routines for the virtual detector system. The
// object 'observe' manages a spherical web of detectors. It contains
// routines to register photons escaping from the model region and to
// save the acquired polarimetric data of each detector to a disk.
//
// STOKES, version 1.0, Nov 2004

#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "observe-stokes-v1.0.h"


// initialize detectors by setting all channels to zero
Monitor::Monitor()
{
  Col = new Column[ScreenNumThetaMax];
}


// register photon by the detector in its flight direction
// and save its wavelength and Stokes vector
void Monitor::RegisterPhoton(Photon p,Model M)
{
  int i,j,ThetaNum,PhiNum;
  double Coords[3],phi,theta;

  // get the number of detectors in the theta and in the phi direction
  PhiNum = M.ScreenNumPhi();
  ThetaNum = M.ScreenNumTheta();

  // compute angular difference between neighboring viewing angles
  PhiStep = 2*pi / PhiNum;
  ThetaStep = pi / ThetaNum;

  // get the photon's flight direction
  Coords[0] = p.GetDirect1();
  Coords[1] = p.GetDirect2();
  Coords[2] = p.GetDirect3();
  SphericalCoords(Coords,phi,theta);

  // if the model space is symmetric with respect to the xy-plane
  // combine the results of the two half-spaces
  if (M.GetPlaneSymmetry())
  {
    if (Coords[2] < 0)
    {
	p.TurnDirection(0,pi);
	Coords[2] = -Coords[2];
    }

    // find the correct detector in the theta direction
    i = int(floor((1-Coords[2])*ThetaNum));

  }

  // else treat all theta directions between 0 and pi
  else
  {

    // find the correct detector in the theta direction of the photon
    i = int(floor(0.5*(1-Coords[2])*ThetaNum));

  }

  // correct i if it is off-limits due to numerical inaccuracies
  if (i<0)
  {
     i=0;
  }
  if (i>=ThetaNum)
  {
     i=ThetaNum-1;
  }

  // adjust the photon's polarization plane to the detector
  p.AdjustFrame(phi,theta);

  // find the correct detector in the phi direction
  j = int(floor(phi/PhiStep));
  if (j<0)
  {
    j=0;
  }
  if (j>=PhiNum)
  {
     j=PhiNum-1;
  }

  // register the photon
  Col[i].Row[j].RegisterPhoton(p,M);

}


// save the data of all detectors to a disk
void Monitor::SaveResults(long int Photons,Model M)
{
  int i,j,n;
  double Comp[3];
  char outputfile[80],dustfile[80];
  CTRegion CTR;
  BLRegion BLR;
  NLRegion NLR;
  scatterer Cloud;
  ofstream out_stream[4];
  time_t rawtime;

  // prepare file names for the output files, there is one file for
  // each Stokes parameter
  for(n=0; n<4; n++)
  {
    M.GetOutputFile(outputfile);
    switch(n)
    {
    case 0:
      strcat(outputfile,"I.sto");
      out_stream[n].open(outputfile,ios::out);
      break;
    case 1:
      strcat(outputfile,"Q.sto");
      out_stream[n].open(outputfile,ios::out);
      break;
    case 2:
      strcat(outputfile,"U.sto");
      out_stream[n].open(outputfile,ios::out);
      break;
    default:
      strcat(outputfile,"V.sto");
      out_stream[n].open(outputfile,ios::out);
      break;
    }
    outputfile[strlen(outputfile)-4]='\0';

    // save all parameters of the model run -

    // ...model run, program version, date and time, dust model
    M.GetDustFile(dustfile);
    dustfile[strlen(dustfile)-4]='\0';
    time(&rawtime);
    out_stream[n] << outputfile << " Stokes version 1.0 "
                  << dustfile << " " << ctime(&rawtime) << "\n";

    // ... emitted spectrum:
    out_stream[n] << Photons << " " << M.PhotonNum() << " ";
    out_stream[n] << M.LambdaMin() << " ";
    out_stream[n] << M.LambdaMax() << " ";
    out_stream[n] << M.LambdaRange() << " ";

    // ...continuum emission region
    CTR = M.GetCTR();
    out_stream[n] << CTR.R << " ";
    out_stream[n] << CTR.a << " ";
    out_stream[n] << CTR.b << " ";
    out_stream[n] << CTR.alpha << " ";
    out_stream[n] << CTR.vr << " ";

    // ...broad line emission region
    BLR = M.GetBLR();
    out_stream[n] << BLR.R << " ";
    out_stream[n] << BLR.a << " ";
    out_stream[n] << BLR.b << " ";
    out_stream[n] << BLR.lambda << " ";
    out_stream[n] << BLR.width << " ";
    out_stream[n] << BLR.rli << " ";
    out_stream[n] << BLR.vr << " ";

    // ...narrow line emission region
    NLR = M.GetNLR();
    out_stream[n] << NLR.a << " ";
    out_stream[n] << NLR.b << " ";
    out_stream[n] << NLR.theta << " ";
    out_stream[n] << NLR.lambda << " ";
    out_stream[n] << NLR.width << " ";
    out_stream[n] << NLR.rli << " ";
    out_stream[n] << NLR.vr << " ";

    // ...web detector
    out_stream[n] << M.ScreenNumPhi() << " ";
    out_stream[n] << M.ScreenNumTheta() << " ";
    out_stream[n] << M.res() << " ";
    out_stream[n] << M.GetPlaneSymmetry() << " ";

    // ...dust properties
    out_stream[n] << M.RadiusMin() << " ";
    out_stream[n] << M.RadiusMax() << " ";
    out_stream[n] << M.GetRadAlpha() << " ";
    M.GetComp(Comp);
    for(j=0; j<=2; j++) out_stream[n] << Comp[j] << " ";
    out_stream[n] << M.GetRadNum() << " ";
    out_stream[n] << M.GetLamNum() << " ";

    // ...computational parameters
    out_stream[n] << M.GetGap() << " ";

    // save data of the scattering regions defined
    out_stream[n] << M.GetScattNum() << " ";
    for(i=0; i<M.GetScattNum(); i++)
    {
      Cloud = M.GetCloud(i);
      out_stream[n] << Cloud.shape << " ";
      switch(Cloud.shape)
      {
        case 0:
          out_stream[n] << Cloud.R << " ";
          out_stream[n] << Cloud.a << " ";
          out_stream[n] << Cloud.b << " ";
          break;
        case 1:
          out_stream[n] << Cloud.R << " ";
          out_stream[n] << Cloud.a << " ";
          out_stream[n] << Cloud.b << " ";
          break;
        case 2:
          out_stream[n] << Cloud.a << " ";
          out_stream[n] << Cloud.b << " ";
          out_stream[n] << 180/pi*Cloud.theta << " ";
          break;
        default:
          for(j=0; j<=2; j++) out_stream[n] << Cloud.pos[j] << " ";
          out_stream[n] << Cloud.a << " ";
          break;
      }
      out_stream[n] << Cloud.dust << " ";
      out_stream[n] << Cloud.Density << " ";
      out_stream[n] << Cloud.vr << " ";
    }

    // save tables of the Stokes parameters -
    out_stream[n] << "\n\n*****\n\n";

    // vertical: viewing direction in theta and phi
    for (i=0; i<M.ScreenNumTheta(); i++)
    {
      for (j=0; j<M.ScreenNumPhi(); j++)
      {
      // horizontal: spectroscopic channels
      for(int itmp=0; itmp < M.res(); itmp++)
        {
          switch(n)
          {
            case 0:
              out_stream[n] << Col[i].Row[j].I[itmp] << " ";
              break;
            case 1:
              out_stream[n] << Col[i].Row[j].Q[itmp] << " ";
              break;
            case 2:
              out_stream[n] << Col[i].Row[j].U[itmp] << " ";
              break;
            default:
              out_stream[n] << Col[i].Row[j].V[itmp] << " ";
              break;
          }
        }
        out_stream[n] << "\n";
      }
      out_stream[n] << "\n";
    }
    out_stream[n].close();
  }
}