;; This buffer is for notes you don't want to save, and for Lisp evaluation.
;; If you want to create a file, visit that file with C-x C-f,
;; then enter the text in that file's own buffer.

// file: seeing.cpp
// ----------------
//
// 

#include "astromath.h"
#include "globalvars.h"
#include "errorc.h"
#include "photometry.h"
#include "coord.h"
#include "seeing.h"


TSeeing::TSeeing()
{
}


TSeeing::~TSeeing()
{
}




// Input parameters:
// D = diameter of entrance pupil (in meters)
// d = distance between entrance pupils (in meters)
// lambda = effective wavelength (in meters)
// VarianceLong = longitudinal variance (in degrees)

double TSeeing::SeeingLong(double VarianceLong, double D, double d, double lambda)
{
 double Varl = VarianceLong; //pow(degrad(sqrt(VarianceLong)),2);
 double r0=pow((2*lambda*lambda*(0.179*pow(D,-1.0/3.0)-0.0968*pow(d,-1.0/3.0))/Varl),3.0/5.0);
 double FWHML=raddeg(0.98*lambda / r0);
 return(FWHML);
}


// Input parameters:
// D = diameter of entrance pupil (in meters)
// d = distance between entrance pupils (in meters)
// lambda = effective wavelength (in meters)
// VarianceTrans = transverse variance (in degrees)

double TSeeing::SeeingTrans(double VarianceTrans, double D, double d, double lambda)
{
 double Vart = VarianceTrans; //pow(degrad(sqrt(VarianceTrans)),2);
 double r0=pow((2*lambda*lambda*(0.179*pow(D,-1.0/3.0)-0.145*pow(d,-1.0/3.0))/Vart),3.0/5.0);
 double FWHMT=raddeg(0.98*lambda / r0);
 return(FWHMT);
}







void TSeeing::CalcSeeing(TStatistics &Statistics,
                         double *FWHMHorL, double *FWHMHorT,
                         double *FWHMVerL, double *FWHMVerT)
{
// D = diameter of entrance pupil (in meters)
// d = distance between entrance pupils (in meters)
// lambda = effective wavelength (in meters)
// VarianceTrans = transverse variance (in degrees)
// PixelHeight = height of image pixels in arcseconds!
// PixelWidth = width of image pixels in arcseconds!

 double D=Statistics.CCDImageProperties.PupilWidth*1e-3;
 double d=Statistics.CCDImageProperties.PupilDistance*1e-3;
 double lambda=Statistics.CCDImageProperties.EffectiveWavelength*1e-9;
// double PixelWidth=Statistics.PixelWidthArcSecs;
// double PixelHeight=Statistics.PixelHeightArcSecs;
 double ZenithDistance=90.0-Statistics.CCDImageProperties.AltAz.Alt();
 double Latitude=atof(GlobalConfiguration.GetCfg("OBSERVATORY","latitude","28.75833333"));; // default La Palma
 

 double SigmaHorX=Statistics.SigmaHorX;
 double SigmaHorY=Statistics.SigmaHorY;
 double SigmaVerX=Statistics.SigmaVerX;
 double SigmaVerY=Statistics.SigmaVerY;
 double MeanXLeft=Statistics.MeanXLeft;
 double MeanXRight=Statistics.MeanXRight;
 double MeanYLeft=Statistics.MeanYLeft;
 double MeanYRight=Statistics.MeanYRight;
 double MeanXUpper=Statistics.MeanXUpper;
 double MeanXBottom=Statistics.MeanXBottom;
 double MeanYUpper=Statistics.MeanYUpper;
 double MeanYBottom=Statistics.MeanYBottom;

 double alpha,beta;
 if (fabs(MeanXRight-MeanXLeft)==0)
   alpha=PI/2;
  else
   alpha=atan(fabs(MeanYRight-MeanYLeft)/fabs(MeanXRight-MeanXLeft));
 if (fabs(MeanYBottom-MeanYUpper)==0)
   beta=PI/2;
  else
   beta=atan(fabs(MeanXBottom-MeanXUpper)/fabs(MeanYBottom-MeanYUpper));

 // this code doesnt' work:
 double SigmaLongHor=SigmaHorX*cos(alpha)+SigmaHorY*sin(alpha);
 double SigmaTransHor=-SigmaHorX*sin(alpha)+SigmaHorY*cos(alpha);
 double SigmaLongVer=SigmaVerY*cos(beta)+SigmaVerX*sin(beta);
 double SigmaTransVer=-SigmaVerY*sin(beta)+SigmaVerX*cos(beta);


 // for the time being, take the old values...
 SigmaLongHor=SigmaHorX;
 SigmaTransHor=SigmaHorY;
 SigmaLongVer=SigmaVerX;
 SigmaTransVer=SigmaVerY;


 *FWHMHorL=SeeingLong(pow(SigmaLongHor,2),D,d,lambda)*3600;
 *FWHMHorT=SeeingTrans(pow(SigmaTransHor,2),D,d,lambda)*3600;
 *FWHMVerL=SeeingLong(pow(SigmaTransVer,2),D,d,lambda)*3600;
 *FWHMVerT=SeeingTrans(pow(SigmaLongVer,2),D,d,lambda)*3600;


 AdjustForZenithAngle(FWHMHorL,FWHMHorT,FWHMVerL,FWHMVerT,  ZenithDistance);
 LogMsg(LOG_Info,"FWHMHorL=%f  FWHMHorT=%f  FWHMVerL=%f  FWHMVerT=%f ",*FWHMHorL,*FWHMHorT,*FWHMVerL,*FWHMVerT);
}




// r0(zenith) = r0(z)  * sec(z)^0.6

void TSeeing::AdjustForZenithAngle(double *FWHMHorL, double *FWHMHorT,
                                   double *FWHMVerL, double *FWHMVerT,
                                   double ZenithDistance) 	// degrees
{
 *FWHMHorL*=pow(1/cos(degrad(ZenithDistance)),0.6);
 *FWHMHorT*=pow(1/cos(degrad(ZenithDistance)),0.6); 
 *FWHMVerL*=pow(1/cos(degrad(ZenithDistance)),0.6);
 *FWHMVerT*=pow(1/cos(degrad(ZenithDistance)),0.6); 
}