#include <stdio.h>
#include <math.h>
#include "SoPo.h"
#include "../General.h"


SoPo::SoPo(float Latitude, float Longitude)
{
    Lat = Latitude  * DEG_TO_RAD;
    Lon = Longitude * DEG_TO_RAD;

    SinLat = sin(Lat);
    CosLat = cos(Lat);
}

void SoPo::Calc_SoPo(Date DD) 
{
 
    
//   float L0,M,e,C,L_true,f,R,GrHrAngle,Obl,RA,Decl,HrAngle; //elev,azimuth;
   float L0,M,C,L_true,GrHrAngle,Obl,RA,Decl,HrAngle; //elev,azimuth;
    
    
//L0=DEG_TO_RAD*fmod(280.46645+36000.76983*DD.T,360);
// M=DEG_TO_RAD*fmod(357.5291+35999.0503*DD.T,360);

L0 = fmod(4.895062287 + 628.331876*DD.T,TWOPI);
M =  fmod(6.240059066 + 628.301865*DD.T,TWOPI);


//e=0.016708617-0.000042037*DD.T;

C = DEG_TO_RAD*((1.9146-0.004847*DD.T)*sin(M)+(0.019993-0.000101*DD.T)*sin(2*M)+0.00029*sin(3*M));

//f=M+C;

// Obl = DEG_TO_RAD*(23 + 26/60.0 + 21.448/3600.0 - (46.815/3600)*DD.T);
Obl = DEG_TO_RAD*(23.4392911 - 0.013004167*DD.T);

GrHrAngle = 280.46061837 + (360*DD.JDx)%360 + 0.98564736629*DD.JDx + 360.98564736629*DD.JD_frac;

GrHrAngle = fmod(GrHrAngle,360.0);

L_true=fmod(C+L0,TWOPI);

// R=1.000001018*(1-e*e)/(1+e*cos(f));

RA=atan2(sin(L_true)*cos(Obl),cos(L_true));

Decl=asin(sin(Obl)*sin(L_true));

HrAngle = DEG_TO_RAD*GrHrAngle+Lon-RA;

// Calculate Elevation
//elev = RAD_TO_DEG * asin(sin(Lat)*sin(Decl)+CosLat*(cos(Decl)*cos(HrAngle)));
r_elev =  asin(SinLat*sin(Decl)+CosLat*(cos(Decl)*cos(HrAngle)));

//===========================0===================
// solar parallax: the difference in position of the Sun as seen from the Earth's centre 
// and a point one Earth radius away, i. e., the angle subtended at the Sun by the Earth's mean radius. 
// Correct for parallax, Horizontal parallax = 8.794" = 4.2635e-5 rad
//r_elev = r_elev - 4.2635e-5 * cos(r_elev);

deg_elev = RAD_TO_DEG * r_elev;

// Azimuth measured eastward from north.
//azimuth = RAD_TO_DEG * (PI+atan2(sin(HrAngle),cos(HrAngle)*SinLat-tan(Decl)*CosLat));
r_azimuth = (PI+atan2(sin(HrAngle),cos(HrAngle)*SinLat-tan(Decl)*CosLat));
deg_azimuth = RAD_TO_DEG * r_azimuth;      
    
// Formula in degrees,  Correct for atmospheric refraction, https://plantpredict.com/algorithm/solar-position/
//volatile   float d_elev = (p0/1010) * (283/(Temp + 273.15) ) * 1.02 / (60*tan(r_elev + (10.3/(deg_elev + 5.11))  )) ;
volatile   float d_elev = (p0/1010) * (283/(Temp + 273.15) ) * 1.02 / (60*tan(r_elev + (DEG_TO_RAD*10.3/(deg_elev + 5.11))  )) ;

 //  elev_corr = r_elev * RAD_TO_DEG;
    deg_elev_corr = deg_elev + d_elev;


}

//========================================================================

void Date::SetDate(int y, int m, int d)
{
    year = y;
    month = m;
    day = d;
    
     if (month<=2) 
        {
            year--; month+=12;
        }
    
// float t0 = (float)year - 2000.0;
//  Delta_T = 62.92 + 0.32217 * t0 + 0.005589 * t0 * t0;
//  Delta_T = 0;
    
}


//========================================================================

void Date::CalcJulian(int h,  int mm, int s)
 {
   
    hour = h;
    minute = mm;
    second = s;
 
// https://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
// https://en.wikipedia.org/wiki/%CE%94T#/media/File:Delta_t.svg

// TT = UT + ?T , TT = terrestrial time (needed), UT = universal time (given by clock)
//    second_deltaT = second + Delta_T; // correct for delta_T till 2050
          
       
    int A,B;

    A=year/100; B=2-A+A/4;
    JD_whole=(int32_t)(365.25*(year+4716))+(int)(30.6001*(month+1))+day+B-1524;
    
     JD_frac=(hour + minute/60.0 + second/3600.0)/24.0 -0.5;
//    JD_frac=(hour + minute/60.0 + second_deltaT/3600.0)/24.0 -0.5;
    T = JD_whole-2451545; 
    
    T = (T+JD_frac)/36525.0;

//    volatile float TT = T;
    
    JDx = JD_whole-2451545;
 }

//========= Create Date and  SoPo instance ================

   //   Date(int y, int m, int d,int h,  int mm, int s)