// Implementation file for the 'model' routines. The object 'model'
// reads in and stores all parameters for the model run defined in
// 'input.par'. Basic parameters and structures of the model space,
// such as emission and scattering regions, are defined here. The
// routines of the object 'model' read in the input data from the file
// 'input.par' and compute basic parameters of the specified model run.
//
// STOKES, version 1.0, Nov 2004


#include "mie-stokes-v1.0.h"
#include "model-stokes-v1.0.h"
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>


// Initialize the Model parameters by defining an empty model
// space. Assume a flat intensity spectrum over the whole possible
// wavelength range of 1000 Ang to 10000 Ang coming from a point-like
// continuum emission region. The (not used) line emission regions are
// set to a cylindrical broad line region and to  a double conical
// narrow line region for the H-alpha line at 6563 Ang. The (not used)
// dust composition reproduces the Galactic extinction curve. In the
// theta direction ten viewing angles are defined, a spectroscopic
// view with 45 channels is declared. A million photons are sampled.
Model::Model()
{
  int i,j;

  PNum = 1000000;
  LamMin = 1000;
  LamMax = 10000;
  CTR.R = 0;
  CTR.a = 0;
  CTR.b = 0;
  CTR.alpha = 2;
  CTR.vr = 0;
  BLR.R = 0.0070;
  BLR.a = 0.0050;
  BLR.b = 0.0050;
  BLR.lambda = 6563;
  BLR.width = 150;
  BLR.rli = 0;
  BLR.vr = 0;
  NLR.a = 10;
  NLR.b = 100;
  NLR.theta = pi/180*25;
  NLR.lambda = 6563;
  NLR.width = 25;
  NLR.rli = 0;
  NLR.vr = 0;
  for(i=0; i<ScattMax; i++)
  {
    Cloud[i].shape = 0;
    for(j=0; j<=2; j++) Cloud[i].pos[j] = 0;
    Cloud[i].R = 0;
    Cloud[i].a = 0;
    Cloud[i].b = 0;
    Cloud[i].theta = 0;
    Cloud[i].vr = 0;
    Cloud[i].dust = false;
    Cloud[i].Density = 0;
  }
  EndScattList = 0;
  RadMin = 0.005;
  RadMax = 0.250;
  RadAlpha = -3.5;
  LamNum = 100;
  RadNum = 100;
  Comp[0] = 0.125;
  Comp[1] = 0.250;
  Comp[2] = 0.625;
  Gap = 1.0e-6;
  SNumPhi = 1;
  SNumTheta = 10;
  planesym = false;
  resolut = 45;
  strcpy(outputfile,"empty-space1");
  strcpy(dustfile,"default.dst");

  // 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 << "                 STOKES\n\n";
  cout << "  Radiative Transfer with Polarization\n\n";
  cout << "              Version 1.0\n";
  cout << "                Nov 2004\n\n";
  cout << "        written by Rene Goosmann\n";
  cout << "    manual: www.stokes-program.info\n";
  cout << "   contact: admin@stokes-program.info\n";
  cout << "========================================\n\n";
}


// read in model parameters from 'input.par', if the initialization
// happens correctly return 'true'.
bool Model::Initialize(Mie& Dust)
{
  bool accept,DustFlag;
  int i,j,FileNum;
  double R,a,b,d;
  char str[200],MStr[4];
  fstream in_stream;

  in_stream.open("input.par",ios::in);
  if (in_stream.fail()!=0)
  {
    cout << "Error - unable to open file 'input.par'!\n";
    cout << "Program aborted...\n\n\n";
    return false;
  }
  else
  {

    // assume that the model does not involve dust
    DustFlag = false;
    
    do
    {
      in_stream >> str;

      // in case of reaching end of file
      if (in_stream.eof())
      {
        strcpy(str,"end");
      }

      // first assume 'str' is not a key word
      accept = false;

      // in case 'str' is a comment go to next line
      if (str[0]=='#')
      {
        in_stream.getline(str,200);
	accept = true;
      }

      // else read in the key word and associated data
      else
      {
        if (strcmp(str,"PhotonNum")==0)
        {
          in_stream >> str;
          PNum = atol(str);
          accept = true;
        }
        if (strcmp(str,"LambdaMin")==0)
        {
          in_stream >> str;
          LamMin = atof(str);
          accept = true;
        }
        if (strcmp(str,"LambdaMax")==0)
        {
          in_stream >> str;
          LamMax = atof(str);
          accept = true;
        }
        if (strcmp(str,"ContSource")==0)
        {
          in_stream >> str;
          CTR.R = atof(str);
          in_stream >> str;
          CTR.a = atof(str);
          in_stream >> str;
          CTR.b = atof(str);
          in_stream >> str;
          CTR.alpha = atof(str);
	  in_stream >> str;
	  CTR.vr = atof(str);
          accept = true;
        }
        if (strcmp(str,"BLSource")==0)
        {
          in_stream >> str;
          BLR.R = atof(str);
          in_stream >> str;
          BLR.a = atof(str);
          in_stream >> str;
          BLR.b = atof(str);
          in_stream >> str;
          BLR.lambda = atof(str);
          in_stream >> str;
          BLR.width = 0.5*atof(str);
          in_stream >> str;
          BLR.rli = 0.01*atof(str);
	  in_stream >> str;
	  BLR.vr = atof(str);
          accept = true;
        }
        if (strcmp(str,"NLSource")==0)
        {
          in_stream >> str;
          NLR.a = atof(str);
          in_stream >> str;
          NLR.b = atof(str);
          in_stream >> str;
          NLR.theta = pi/180*atof(str);
          in_stream >> str;
          NLR.lambda = atof(str);
          in_stream >> str;
          NLR.width = 0.5*atof(str);
          in_stream >> str;
          NLR.rli = 0.01*atof(str);
	  in_stream >> str;
	  NLR.vr = atof(str);
          accept = true;
        }
        if (strcmp(str,"Torus")==0)
        {
          i = EndScattList;
          Cloud[i].shape = 0;
          in_stream >> str;
          Cloud[i].R = atof(str);
          in_stream >> str;
          Cloud[i].a = atof(str);
          in_stream >> str;
          Cloud[i].b = atof(str);
          in_stream >> str;
          Cloud[i].dust = (strcmp(str,"dust")==0);
          if (Cloud[i].dust) DustFlag = true;
          in_stream >> str;
          Cloud[i].Density = atof(str);
          in_stream >> str;
          Cloud[i].vr = atof(str);
          EndScattList++;
          accept = true;
        }
        if (strcmp(str,"Cylinder")==0)
        {
          i = EndScattList;
          Cloud[i].shape = 1;
          in_stream >> str;
          Cloud[i].R = atof(str);
          in_stream >> str;
          Cloud[i].a = atof(str);
          in_stream >> str;
          Cloud[i].b = atof(str);
          in_stream >> str;
          Cloud[i].dust = (strcmp(str,"dust")==0);
          if (Cloud[i].dust) DustFlag = true;
          in_stream >> str;
          Cloud[i].Density = atof(str);
          in_stream >> str;
          Cloud[i].vr = atof(str);
          EndScattList++;
          accept = true;
        }
        if (strcmp(str,"Double-cone")==0)
        {
          i = EndScattList;
          Cloud[i].shape = 2;
          in_stream >> str;
          Cloud[i].a = atof(str);
          in_stream >> str;
          Cloud[i].b = atof(str);
          in_stream >> str;
          Cloud[i].theta = pi/180*atof(str);
          in_stream >> str;
          Cloud[i].dust = (strcmp(str,"dust")==0);
          if (Cloud[i].dust) DustFlag = true;
          in_stream >> str;
          Cloud[i].Density = atof(str);
          in_stream >> str;
          Cloud[i].vr = atof(str);
          EndScattList++;
          accept = true;
        }
        if (strcmp(str,"Sphere")==0)
        {
          i = EndScattList;
          Cloud[i].shape = 3;
          for (j=0; j<=2; j++)
          {
            in_stream >> str;
            Cloud[i].pos[j] = atof(str);
          }
          in_stream >> str;
          Cloud[i].a = atof(str);
          in_stream >> str;
          Cloud[i].dust = (strcmp(str,"dust")==0);
          if (Cloud[i].dust) DustFlag = true;
          in_stream >> str;
          Cloud[i].Density = atof(str);
          in_stream >> str;
          Cloud[i].vr = atof(str);
          EndScattList++;
          accept = true;
        }
        if (strcmp(str,"GrainRadMin")==0)
        {
          in_stream >> str;
          RadMin = atof(str);
          accept = true;
        }
	if (strcmp(str,"GrainRadMax")==0)
        {
	  in_stream >> str;
          RadMax = atof(str);
	  accept = true;
        }
	if (strcmp(str,"GrainSizeInd")==0)
        {
          in_stream >> str;
          RadAlpha = atof(str);
          accept = true;
        }
	if (strcmp(str,"GrainRadNum")==0)
        {
          in_stream >> str;
          RadNum = atoi(str);
          accept = true;
        }
	if (strcmp(str,"GrainLambdaNum")==0)
        {
          in_stream >> str;
          LamNum = atoi(str);
          accept = true;
        }
        if (strcmp(str,"DustComp")==0)
        {
          for (i=0; i<=2; i++)
          {
            in_stream >> str;
            Comp[i] = atof(str)/100;
          }
          accept = true;
        }
        if (strcmp(str,"Gap")==0)
        {
          in_stream >> str;
          Gap = atof(str);
          accept = true;
        }
        if (strcmp(str,"PhiViewAng")==0)
        {
          in_stream >> str;
          SNumPhi = atoi(str);
          accept = true;
        }
        if (strcmp(str,"ThetaViewAng")==0)
        {
          in_stream >> str;
          SNumTheta = atoi(str);
          accept = true;
        }
        if (strcmp(str,"PlaneSym")==0)
        {
          in_stream >> str;
          planesym = (strcmp(str,"yes")==0);
          accept = true;
        }
        if (strcmp(str,"SpectRes")==0)
        {
          in_stream >> str;
          resolut = atoi(str);
          accept = true;
        }
	if (strcmp(str,"OutputFile")==0)
	{
	  in_stream >> outputfile;
	  in_stream >> str;
	  FileNum = atoi(str);
	  sprintf(MStr,"%d",FileNum);
	  strcat(outputfile,MStr);
	  accept = true;
	}
	if (strcmp(str,"DustModel")==0)
	{
	  in_stream >> dustfile;
	  strcat(dustfile,".dst");
	  accept = true;
	}
      }
    }
    while((accept)&&(strcmp(str,"end")!=0));
    in_stream.close();

    // If data read in correctly compute some global parameters
    if (strcmp(str,"end")==0)
    {

      // Determine the maximum expansion among all scatterers and their
      // maximum radial velocity. This provides an upper limit of the
      // overall size and dynamics of the model space as a whole.
      MaxExp = 0;
      VMax = 0;
      for(i=0; i<EndScattList; i++)
      {
	switch(Cloud[i].shape)
	{
	  case 0:
	    a = Cloud[i].a;
	    b = Cloud[i].b;
	    R = Cloud[i].R;
	    d = 2*(sqrt((R+b)*(R+b)+a*a))+Gap;
	    if (d > MaxExp) MaxExp = d;
	    break;
          case 1:
	    a = Cloud[i].a;
	    b = Cloud[i].b;
  	    R = Cloud[i].R;
  	    d = 2*(sqrt((R+b)*(R+b)+a*a))+Gap;
  	    if (d > MaxExp) MaxExp = d;
  	    break;
          case 2:  d = 2*Cloud[i].b+Gap;
    	    if (d > MaxExp) MaxExp = d;
  	    break;
          default: d = 2*(sqrt(SProd(Cloud[i].pos,Cloud[i].pos))+Cloud[i].a)
                       +Gap;
    	    if (d > MaxExp) MaxExp = d;
        }
        if (Cloud[i].vr>VMax)
        {
  	  VMax = Cloud[i].vr;
        }
      }

      // Determine parts of the sampling integrals for the line photons
      // in order to shorten the computations in the main loop
      BLR.low = 0.5*log(pow(LamMin-BLR.lambda,2)+pow(BLR.width,2))+
                BLR.lambda/BLR.width*atang(LamMin-BLR.lambda,BLR.width);
      BLR.norm = 0.5*log(pow(LamMax-BLR.lambda,2)+pow(BLR.width,2))+
                 BLR.lambda/BLR.width*atang(LamMax-BLR.lambda,BLR.width)-
                 BLR.low;
      NLR.low = 0.5*log(pow(LamMin-NLR.lambda,2)+pow(NLR.width,2))+
                NLR.lambda/NLR.width*atang(LamMin-NLR.lambda,NLR.width);
      NLR.norm = 0.5*log(pow(LamMax-NLR.lambda,2)+pow(NLR.width,2))+
                 NLR.lambda/NLR.width*atang(LamMax-NLR.lambda,NLR.width)-
                 NLR.low;

      cout << "Input parameters properly read in\n";
      cout << "Output file name set to '" << outputfile << "'\n";

      // if the model run involves dust scattering compute or read the
      // Mie data
      if (DustFlag) Dust.InitMie(RadMin,RadMax,RadNum,LamNum,
                                 Comp,RadAlpha,dustfile);

      return true;
    }

    else
    {
      cout << "Error - '" << str << "' is not a key word!\n";
      cout << "Program aborted...\n\n\n";
      return false;
    }
  }
}


// return the number of detectors in the phi-direction
int Model::ScreenNumPhi()
{
  return SNumPhi;
}


// return the number of detectors in the theta-direction
int Model::ScreenNumTheta()
{
  return SNumTheta;
}


// return the number of spectroscopic channels
int Model::res()
{
  return resolut;
}


// Return the number of photons to sample
long int Model::PhotonNum()
{
  return PNum;
}


// return the number of scattering regions defined
int Model::GetScattNum()
{
  return EndScattList;
}


// return parameters of scattering region with index i
scatterer Model::GetCloud(int i)
{
  return Cloud[i];
}


// return the maximum expansion of all scattering regions
double Model::GetMaxExp()
{
  return MaxExp;
}


// return the maximum radial velocity of all scattering regions
double Model::GetVMax()
{
  return VMax;
}


// return the minimum photon wavelength defined
double Model::LambdaMin()
{
  return LamMin;
}


// return the maximum photon wavelength defined
double Model::LambdaMax()
{
  return LamMax;
}


// return the photon wavelength range covered
double Model::LambdaRange()
{
  return LamMax - LamMin;
}


// set the output file name to name
void Model::GetOutputFile(char name[80])
{
  strcpy(name,outputfile);
}


// set the dust file name to name
void Model::GetDustFile(char name[80])
{
  strcpy(name,dustfile);
}


// return the maximum grain radius defined
double Model::RadiusMax()
{
  return RadMax;
}


// return the minimum grain radius defined
double Model::RadiusMin()
{
  return RadMin;
}


// return the index of the grain size distribution
double Model::GetRadAlpha()
{
  return RadAlpha;
}


// return the number of grain radii considered for calculating
// the Mie data
int Model::GetRadNum()
{
  return RadNum;
}


// return the number of wavelengths considered for calculating
// the Mie data
int Model::GetLamNum()
{
  return LamNum;
}


// set the percentages of all dust components to 'Components'
void Model::GetComp(double Components[SubNum])
{
  int s;

  for(s=0; s<SubNum; s++) Components[s] = Comp[s];
}


// return the minimum gab between scattering clouds
double Model::GetGap()
{
  return Gap;
}


// return 'true' if the model space is defined symmetrically with
// respect to the xy-plane and 'false' otherwise
bool Model::GetPlaneSymmetry()
{
  return planesym;
}


// return the parameters of the continuum emission region
CTRegion Model::GetCTR()
{
  return CTR;
}


// return the parameters of the broad line emission region
BLRegion Model::GetBLR()
{
  return BLR;
}


// return the parameters of the narrow line emission region
NLRegion Model::GetNLR()
{
  return NLR;
}