// Main part of the program ANALYZE designed to convert STOKES results
// to a convenient format. The program computes the total photon
// number, the degree of polarization, and its position angle. It also
// produces a new set of combined STOKES results files in case one
// wants to expand the simulation in the future by performing further
// runs.
//
// ANALYZE, version 1.0, Nov 2004


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


// definition of pi
const double pi = 3.1415926535898;

// maximum number of model parameters stored in the input files
const int DataMax = 200;

int m,n,InputSets,DataPtr;
double Data[DataMax];
char Inputfile[80],Outputfile[80],MStr[4],infile[80],outfile[80],
     resfile[80],title[10][80];
bool flag;
fstream par_stream,in_stream[4];
fstream out_stream[4],res_stream[3];
Monitor M;


// If one wants to combine STOKES results obtained on different CPUs
// the program checks if the parameters of all output files are
// consistent; if this is not the case an error message is displayed
void Warning(int m,int n,int k)
{ 
  char ch;

  cout << "Warning - data files are not compatible!\n";
  cout << "File ..." << m;
  switch (n)
  {
  case 0:
    cout << "I.sto";
    break;
  case 1:
    cout << "Q.sto";
    break;
  case 2:
    cout << "U.sto";
    break;
  default:
    cout << "V.sto";
    break;
  }
  if (k==0)
  {
    cout << " has differing model parameters.\n";
  }
  else
  {
    cout << " was not created by version 1.0 of STOKES.\n";
  }
}


// read the data file name from 'input.par' and initialize the
// parameter array by setting all values to zero; return 'true' if no
// error occurs, and 'false' if 'input.par' cannot be opened
bool Initialize(char Inputfile[80],char Outputfile[80],
               double Data[DataMax],int DataPtr)
{
  int i;
  char str[200];

  // set output format for numbers
  cout.setf(ios::fixed);
  cout.setf(ios::showpoint);
  cout.precision(5);

  // welcome screen
  cout << "\n\n";
  cout << "========================================\n\n";
  cout << "                ANALYZE\n\n";
  cout << "   Data Processing of STOKES Results\n\n";
  cout << "              Version 1.0\n";
  cout << "                Nov 2004\n\n";
  cout << "========================================\n\n";
  par_stream.open("input.par",ios::in);
  if (par_stream.fail()!=0)
  {
    return false;
  }
  else
  {
    Inputfile[0]='\0';
    while (!par_stream.eof())
    {
      par_stream >> str;
      if (str[0]=='#')
      {
        par_stream.getline(str,200);
      }
      else
      {
        if (strcmp(str,"OutputFile")==0)
        {
          par_stream >> Inputfile;
          break;
        }
      }
    }
    if (Inputfile[0]=='\0') strcpy(Inputfile,"empty-space");
    strcpy(Outputfile,Inputfile);
    for(i=0; i<DataMax; i++)
    {
      Data[i] = 0;
    }
    DataPtr = 0;
    return true;
  }
}


// the program reads data files created by STOKES and produces
// new output files; here the required suffix for each case is added
// to the file name
void AttachSuffix(char filename[80],int n,int type)
{
  n = 4*type+n;
  switch(n)
  {
  case 0:
    strcat(filename,"I.sto");
    break;
  case 1:
    strcat(filename,"Q.sto");
    break;
  case 2:
    strcat(filename,"U.sto");
    break;
  case 3:
    strcat(filename,"V.sto");
    break;
  case 4:
    strcat(filename,"TF.res");
    break;
  case 5:
    strcat(filename,"PO.res");
    break;
  case 6:
    strcat(filename,"PA.res");
    break;
  default:
    strcat(filename,"PF.res");
    break;
  }
}


// read in the title elements of the first set of files
void GetTitle(fstream& in_stream,char title[10][80],int m,int n)
{
  int i;
    
  for(i=0; i<10; i++)
  {
    in_stream >> title[i];
  }
  title[0][strlen(title[0])-2] = '\0';
  if (strcmp(title[3],"1.0")!=0) Warning(m,n,1);
}


// read in the parameters of the first set of files, compare the
// parameters of all other sets with the first one; if there is an
// inconsistency print a warning message
void GetData(fstream& in_stream,double Data[DataMax],int& DataPtr,
             int m,int n)
{
  int i;
  char input[80];
  bool flag,diffparam;

  i = 0;
  flag = false;
  diffparam = false;
  while(!flag)
  {
    in_stream >> input;

    // read until end of parameter section
    if (strcmp(input,"*****")==0)
    {
      flag = true;
    }
    
    else
    {

      // add up total number of photons sampled/to sample
      if (i<2)
      {
        if (n==0)
        {
          Data[i] = Data[i]+atof(input);
        }
      }

      // read in data of the first set of files, compare
      // the parameters of all other files to this set
      else
      {
        if ((m==1)&&(n==0))
        {
          Data[i] = atof(input);
        }
        else
        {
          if((atof(input)!=Data[i])&&(!diffparam))
          {
            Warning(m,n,0);
	    diffparam = true;
          }
        }
      }
      i++;
    }
  }
  DataPtr = i;
}


// add up the Stokes parameters of all files
void AddValues(fstream& in_stream,double Data[],int m,int n,Monitor M)
{
  int i,j,l;
  double value;

  // loop over viewing angle theta
  for(i=0; i<Data[25]; i++)
  {

    // loop over viewing angle phi
    for(j=0; j<Data[24]; j++)
    {

      // loop over spectral channels
      for(l=0; l<Data[26]; l++)
      {

        // read value
        in_stream >> value;

	// add value to the correct data set
        M.AddTo(i,j,n,l,value);

      }
    }
  }
}


// write all parameters to a set of combined STOKES data files, for
// future expansion of the modeling, and to a set of output files
// containing polarization data
void SaveParameters(char outfile[80],char resfile[80],double Data[],
                    int DataPtr,int InputSets,char title[10][80],
                    int n,Monitor M)
{
  int i,j;
  fstream out_stream,res_stream;
  time_t rawtime;

  // create combined data files and output files
  out_stream.open(outfile,ios::out);
  res_stream.open(resfile,ios::out);

  // save title elements for combined data files
  out_stream << title[0];
  switch(n)
  {
    case 0:
      out_stream << "CompI ";
      break;
    case 1: 
      out_stream << "CompQ ";
      break;
    case 2: 
      out_stream << "CompU ";
      break;
    default:
      out_stream << "CompV ";
      break;
  }
  for(i=1; i<10; i++)
  {
    out_stream << title[i] << " ";
  }
  out_stream << "\n\n";
  
  // save parameters for the combined data files
  for(i=0; i<DataPtr; i++)
  {
    out_stream << Data[i] << " ";
  }
  out_stream << "\n\n*****\n\n";

  // save parameters for the output files:

  // general model data:
  switch(n)
  {
    case 0:
      res_stream << "Total flux\n\n";
      break;
    case 1: 
      res_stream << "Percentage of polarization\n\n";
      break;
    case 2: 
      res_stream << "Polarization position angle\n\n";
      break;
    default:
      res_stream << "Polarized flux\n\n";
      break;
  }
  res_stream << "Model: " << title[0] << "\n" << "      ";
  res_stream << "Data-sets " << InputSets << "\n" << "      ";
  res_stream << "Stokes-version " << title[3] << "\n" << "      ";
  res_stream << "Analyze-version 1.0" << "\n" << "      ";
  time(&rawtime);
  res_stream << "Processed " << ctime(&rawtime) << "\n";
  res_stream << "Photons: " << Data[0] << " / " << Data [1] << "\n\n";
  res_stream << "Lambda: ";
  res_stream << "Minimum " << Data[2] << " Ang\n" << "      ";
  res_stream << "Maximum " << Data[3] << " Ang\n" << "      ";
  res_stream << "Range " << Data[4] << " Ang\n\n";
  res_stream << "Continuum: ";
  res_stream << "R= " << Data[5] << " pc\n      ";
  res_stream << "a= " << Data[6] << " pc\n      ";
  res_stream << "b= " << Data[7] << " pc\n      ";
  res_stream << "Spectral-index " << Data[8] << "\n      ";
  res_stream << "Velocity= " << Data[9] << " km/s\n\n";
  res_stream << "BLR: ";
  res_stream << "R= " << Data[10] << " pc\n      ";
  res_stream << "a= " << Data[11] << " pc\n      ";
  res_stream << "b= " << Data[12] << " pc\n      ";
  res_stream << "Lambda= " << Data[13] << " Ang\n      ";
  res_stream << "Width= " << Data[14]*2 << " Ang\n      ";
  res_stream << "Rel.-strength= " << Data[15]*100 << " %" << "\n      ";
  res_stream << "Velocity= " << Data[16] << " km/s\n\n";
  res_stream << "NLR: ";
  res_stream << "a= " << Data[17] << " pc\n      ";
  res_stream << "b= " << Data[18] << " pc\n      ";
  res_stream << "Theta= " << Data[19]*180/pi << " deg\n      ";
  res_stream << "Lambda= " << Data[20] << " Ang\n      ";
  res_stream << "Width= " << Data[21]*2 << " Ang\n      ";
  res_stream << "Rel.-strength= " << Data[22]*100 << " %" << "\n      ";
  res_stream << "Velocity= " << Data[23] << " km/s\n\n";
  res_stream << "Detectors: ";
  res_stream << "Phi " << Data[24] << "\n" << "      ";
  res_stream << "Theta " << Data[25] << "\n" << "      ";
  res_stream << "Resolution " << Data[26] << " channel(s)\n\n";
  res_stream << "Plane-symmetry: ";
  if (Data[27]==1)
  {
      res_stream << "yes\n\n";
  }
  else
  {
      res_stream << "no\n\n";
  }

  // dust data:
  res_stream << "Dust: " << title[4] << "\n" << "      ";
  res_stream << "GrainRadmin " << Data[28] << " micrometer\n" << "      ";
  res_stream << "GrainRadmax " << Data[29] << " micrometer\n" << "      ";
  res_stream << "GrainSizeIndex " << Data[30] << "\n" << "      ";
  res_stream << "Para-Graphite " << Data[31]*100 << " %\n" << "      ";
  res_stream << "Ortho-Graphite " << Data[32]*100 << " %\n" << "      ";
  res_stream << "Silicate " << Data[33]*100 << " %\n" << "      ";
  res_stream << "bins/radii " << Data[34] << "\n" << "      ";
  res_stream << "bins/wavelength " << Data[35] << "\n\n";
  
  // minimum distance between clouds
  res_stream << "Cloud-gap: " << Data[36] << " pc\n\n";

  // data for the scattering regions
  res_stream << "Scattering-region(s): ";
  res_stream << Data[37] << "\n\n";
  j = 38;
  for(i=0; i<Data[37]; i++)
  {
    res_stream << i+1 << ") ";
    switch(int(floor(Data[j])))
    {
      case 0:
        res_stream << "Torus:" << "\n";
        res_stream << "      " << "R= " << Data[j+1] << " pc\n";
        res_stream << "      " << "a= " << Data[j+2] << " pc\n";
        res_stream << "      " << "b= " << Data[j+3] << " pc\n";
        j=j+4;        
        break;
      case 1:
        res_stream << "Cylinder:" << "\n";
        res_stream << "      " << "R= " << Data[j+1] << " pc\n";
        res_stream << "      " << "a= " << Data[j+2] << " pc\n";
        res_stream << "      " << "b= " << Data[j+3] << " pc\n";
        j=j+4;
        break;
      case 2:
        res_stream << "Double-cone:" << "\n";
        res_stream << "      " << "a= " << Data[j+1] << " pc\n";
        res_stream << "      " << "b= " << Data[j+2] << " pc\n";
        res_stream << "      " << "Theta= " << Data[j+3] << " deg\n";
        j=j+4;
        break;
      default:
        res_stream << "Sphere:" << "\n";
        res_stream << "      " << "(x,y,z)= ( ";
        res_stream << Data[j+1] << " pc, ";
        res_stream << Data[j+2] << " pc, ";
        res_stream << Data[j+3] << " pc )\n";
        res_stream << "      " << "r= " << Data[j+4] << " pc\n";
        j=j+5;
        break;
    }
    res_stream << "      ";
    if(Data[j]==1)
    {
      res_stream << "Dust-density= ";
    }
    else
    {
      res_stream << "Electron-density= ";
    }
    res_stream << Data[j+1] << " /ccm\n"; 
    res_stream << "      ";
    res_stream << "Velocity= " << Data[j+2] << " km/s\n\n";
    j=j+3;
  }

  // close files
  out_stream.close();
  res_stream.close();

}


// main program - this routine reads in all STOKES data files that can
// be found in the current directory, combines the data, and produces
// convenient output files of the polarization results
int main()
{

  // if file 'input.par' found
  if (Initialize(Inputfile,Outputfile,Data,DataPtr))
  {
    strcpy(infile,Inputfile);
    flag = false;
    m = 0;
  
    // search for successive data files
    do
    {
  
      // construct data file name
      m++;
      sprintf(MStr,"%d",m);
  
      // read the data files for the four Stokes parameters
      for(n=0; n<4; n++)
      {
        infile[strlen(Inputfile)] = '\0';
        strcat(infile,MStr);
        AttachSuffix(infile,n,0);
        in_stream[n].open(infile,ios::in);
        if (in_stream[n].fail()!=0)              
        {
          flag = true;
          InputSets = m-1;
          break;
        }
  
        // read data, add it up, and store it in the object observe 
        else
        {
          GetTitle(in_stream[n],title,m,n);
          GetData(in_stream[n],Data,DataPtr,m,n);
          AddValues(in_stream[n],Data,m,n,M);
        }
  
        in_stream[n].close();
      }
      if (!flag) cout << "Read input file set " << Inputfile << " " 
                      << m << "\n";
    }
    while (!flag);
  
    // print error message if first data set is not found
    if (InputSets==0)
    {
      cout << "Error: file set " << Inputfile << "1 not found -\n";
      cout << "Program aborted.\n\n\n";
    }
  
    // else process data files
    else
    {
  
      // construct file names for the combined data files and for the
      // output files
      strcpy(outfile,Outputfile);
      strcpy(resfile,Outputfile);
      for(n=0; n<4; n++)
      {
        outfile[strlen(Outputfile)] = '\0';
        AttachSuffix(outfile,n,0);
        resfile[strlen(Outputfile)] = '\0';
        AttachSuffix(resfile,n,1);
  
        // save parameters and all data
        SaveParameters(outfile,resfile,Data,DataPtr,InputSets,title,n,M);
        M.SaveResults(outfile,resfile,n,Data);
  
      }
      cout << "Created output files\n";
      cout << "Program finished.\n\n\n";
    }
  }

  // file 'input.par' not found
  else
  {
    cout << "Error: unable to open file 'input.par'\n";
    cout << "Program aborted.\n\n\n";
  }

  return 0;
}