// Implementation of the routines of the object 'observe'. This object
// is fairly similar to the one defined in STOKES and manages a web of
// detectors. The program ANALYZE will store the sum of all data files
// in these detectors and finally save their contents into a new set
// of STOKES files and in convenient output files for the polarization
// data.
//
// ANALYZE, version 1.0, Nov 2004


#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include "observe-analyze-v1.0.h"


// definition of pi
const double pi = 3.1415926535898;


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


// add 'value' to the channel l of the detector given by the theta and
// phi position i and j; n denotes the number of the Stokes parameter
void Monitor::AddTo(int i,int j,int n,int l,double value)
{
  Col[i].Row[j].AddTo(n,l,value);
}


// save the data of all detectors to new combined STOKES files and to
// convenient output files
void Monitor::SaveResults(char outfile[20],char resfile[20],
                          int n,double Data[])
{
  int i,j,l;
  double k,dk,N,p,theta;
  double I[ResMax],Q[ResMax],U[ResMax],V[ResMax];
  ofstream out_stream,res_stream;

  // open the files
  out_stream.open(outfile,ios::app);
  res_stream.open(resfile,ios::app);

  // loop for the theta direction
  for (i=0; i<Data[25]; i++)
  {

    // print a line of wavelengths of the channels
    if ((i==0)||(Data[24]>1))
    {
      res_stream << "***** ";
      dk = Data[4]/Data[26];
      for(k=Data[2]+dk/2; k<Data[3]; k=k+dk)
      {
        res_stream << k << " ";
      }
      res_stream << "\n";
    }

    // loop for the phi direction
    for (j=0; j<Data[24]; j++)
    {

      // print cos(theta) for a model space symmetric to the xy-plane
      if (Data[27]==1)
      {
         res_stream << "cosT" << 1 - (i+0.5)/Data[25];
      }

      // print cos(theta) for a model space not symmetric to the
      // xy-plane
      else
      {
         res_stream << "cosT" << 1 - 2*(i+0.5)/Data[25];
      }

      // if there are more than 1 azimthal viewing angles print 
      // the value of phi
      if (Data[24]>1)
      {
        res_stream << "/P" << (j+0.5)*360.0/Data[24] << " ";
      }
	else
      {
        res_stream << " ";
      }

      // loop for the spectroscopic channels
      for (l=0; l<Data[26]; l++)
      {
        I[l] = Col[i].Row[j].I[l];
        Q[l] = Col[i].Row[j].Q[l];
        U[l] = Col[i].Row[j].U[l];
        V[l] = Col[i].Row[j].V[l];

	// compute the polarization p
        if (I[l]!=0)
        {
          p = sqrt(Q[l]*Q[l]+U[l]*U[l]+V[l]*V[l])/I[l];
        }
        else
        {
          p = 0;
        }

        // compute the polarization position angle
        if ((U[l]==0)&&(Q[l]==0))
        {
          theta = 0;
        }
        else
        {
          theta = 90/pi*atan2(U[l],Q[l]);
        }

	// compute a normalization constant for the fluxes
        N = (Data[0]*(Data[2]+(l+0.5)*Data[4]/Data[26]))/10000;

        switch(n)
        {

          // save combined I-parameter and normalized total flux
          case 0:
            out_stream << I[l] << " ";
            res_stream << I[l]/N << " ";
            break;

	  // save combined Q-parameter and polarization
          case 1:
            out_stream << Q[l] << " ";
            res_stream << p << " ";
            break;

	  // save combined U-parameter and position angle
          case 2:
            out_stream << U[l] << " ";
            if (theta < -45) theta = theta + 180;
            res_stream << theta << " ";
            break;

	  // save combined V parameter and normalized polarized flux
          default:
            out_stream << V[l] << " ";
            res_stream << p*I[l]/N << " ";
            break;

        }

      }

      out_stream << "\n";
      res_stream << "\n";
    }

  }

  out_stream.close();
  res_stream.close();
}