//#include <iostream>
#include <inttypes.h>
//typedef char byte;
//#define float double
using namespace std;
#define PROGMEM
#define PI M_PI
#define DEG_TO_RAD (PI/180)
#define pgm_read_float(o) (*(o))

#include <math.h>
//#include "orbitalElements.h"
#include <avr/pgmspace.h>
//http://www.stjarnhimlen.se/comp/ppcomp.html
//sun orbital elements
const PROGMEM  float soe[]  = { 0,0,282.9404,4.70935e-5,1,
                                0.016709,-1.151e-9,356.047,0.9856002585
                              };
#define SOE_N (pgm_read_float(soe))
#define SOE_I (pgm_read_float(soe+1))
#define SOE_W(day) (pgm_read_float(soe+2)+(day)*(pgm_read_float(soe+3)))
#define SOE_A (pgm_read_float(soe+4))
#define SOE_E(day) (pgm_read_float(soe+5)+(day)*(pgm_read_float(soe+6)))
#define SOE_M(day) (pgm_read_float(soe+7)+(day)*(pgm_read_float(soe+8)))

//moon orbital elements
const PROGMEM  float moe[]  = { 125.1228, -0.0529538083, 5.1454, 318.0634, 0.1643573223,
                                60.266,0.0549,115.3654, 13.0649929509
                              };
#define MOE_N(day) (pgm_read_float(moe)+(day)*(pgm_read_float(moe+1)))
#define MOE_I (pgm_read_float(moe+2))
#define MOE_W(day) (pgm_read_float(moe+3)+(day)*(pgm_read_float(moe+4)))
#define MOE_A (pgm_read_float(moe+5))
#define MOE_E (pgm_read_float(moe+6))
#define MOE_M(day) (pgm_read_float(moe+7)+(day)*(pgm_read_float(moe+8)))

float sunRise, sunSet,
      moonRise, moonSet;
double dayNumber(int year, byte m, byte D, double UT)
{
    return (double)(367 * (long)year - 7 * (year + (m + 9) / 12) / 4 + 275 * m / 9 + D - 730530) + UT / 24;
}

void sun(int year, byte m, byte D, double UT)
{
    float h = -0.8 * DEG_TO_RAD;
    float d = dayNumber(year, m, D, UT);
    float M = SOE_M(d);
    while (M >= 360)M -= 360;
    M *= DEG_TO_RAD;
    float e = SOE_E(d);
    float E = M + e * sin(M) * (1.0 + e * cos(M));
    float xv = cos(E) - e;
    float yv = sqrt(1.0 - e * e) * sin(E);
    float v = atan2(yv, xv);
    float r = sqrt(xv * xv + yv * yv);
    float lonsun = ( v + SOE_W(d) * DEG_TO_RAD );
    float xs = r * cos(lonsun);
    float ys = r * sin(lonsun);
    float xe = xs;
    float ecl = (23.4393 - 3.563e-7 * d) * DEG_TO_RAD;
    float ye = ys * cos(ecl);
    float ze = ys * sin(ecl);
    float sunRA = atan2(ye, xe);
    float sunDec = atan2(ze, sqrt(xe * xe + ye * ye));
    float longit = 20 * DEG_TO_RAD; //wspolrzedne geograficzne +20E,+50N
    float lat = 50 * DEG_TO_RAD;
    float UT_Sun_in_south = (sunRA - (PI + M + SOE_W(d) * DEG_TO_RAD) - longit) / (15 * DEG_TO_RAD);
    while (UT_Sun_in_south < 0) UT_Sun_in_south += 24;
    float clha = (sin(h) - sin(lat) * sin(sunDec)) / (cos(lat) * cos(sunDec));
    if (abs(clha) <= 1) clha = acos(clha);
    else
    {
        /*polarny dzien lub noc, lub Slonce nie schodzi ponizej h*/
    }
    sunRise = UT_Sun_in_south - clha * 12 / PI;
    sunSet = UT_Sun_in_south + clha * 12 / PI;
}
void moon(int year, byte m, byte D, double UT)
{
    float h = 0 * DEG_TO_RAD;
    float d = dayNumber(year, m, D, UT);
    float M = MOE_M(d);
    while (M >= 360)M -= 360;
    while (M < 0) M += 360;
    M *= DEG_TO_RAD;
    float e = MOE_E;
    float E = M + e * sin(M) * (1.0 + e * cos(M));
    float a=MOE_A;
    float xv = a*(cos(E) - e);
    float yv = a*sqrt(1.0 - e * e) * sin(E);
    float v = atan2(yv, xv);
    float r = sqrt(xv * xv + yv * yv);
    float N=MOE_N(d);
    while (N>=360) N-=360;
    while (N<0) N+=360;
    N *= DEG_TO_RAD;
    float w=MOE_W(d);
    while (w>=360) w-=360;
    while (w<0) w+=360;
    w *= DEG_TO_RAD;
    float i=MOE_I*DEG_TO_RAD;
    float xh=r*(cos(N)*cos(v+w)-sin(N)*sin(v+w)*cos(i));
    float yh=r*(sin(N)*cos(v+w)+cos(N)*sin(v+w)*cos(i));
    float zh=r*(sin(v+w)*sin(i));
    float lonecl=atan2(yh,xh);
    float latecl=atan2(zh,sqrt(xh*xh+yh*yh));
    float Ls=(SOE_M(d)-SOE_W(d))*DEG_TO_RAD;
    float Lm=M+N+w;
    lonecl+=(-1.274*sin(M-2*(Lm-Ls))+.658*sin(2*(Lm-Ls))-0.186*sin(SOE_M(d)*DEG_TO_RAD))*DEG_TO_RAD;
    latecl+=(-0.173*sin(M+w-2*(Lm-Ls))+0.055*sin(w+2*(Lm-Ls)))*DEG_TO_RAD;//tutaj przy 0.055 troche poskracalem, bo dalo sie :)
    r-=(0.58*cos(M-2*(Lm-Ls))+0.46*cos(2*(Lm-Ls)));
    xh=r*cos(lonecl)*cos(latecl);
    yh=r*sin(lonecl)*cos(latecl);
    zh=r*sin(latecl);
    float ecl= (23.4393 - 3.563E-7 * d)*DEG_TO_RAD;
    float xe=xh;
    float ye=yh*cos(ecl)-zh*sin(ecl);
    float ze=yh*sin(ecl)+zh*cos(ecl);
    float moonRA=atan2(ye,xe);
    float moonDec = atan2(ze, sqrt(xe*xe+ye*ye));
    float longit = 20 * DEG_TO_RAD; //wspolrzedne geograficzne +20E,+50N
    float lat = 50 * DEG_TO_RAD;
// cout << "RA "<< moonRA*12/PI << ":Dec " << (moonDec/(DEG_TO_RAD) )<<" r:" << r<<endl;
    float UT_Moon_in_south = (moonRA - ( PI+(SOE_M(d)+SOE_W(d))*DEG_TO_RAD) - longit) / (15 * DEG_TO_RAD);
    while (UT_Moon_in_south < 0) UT_Moon_in_south += 24;
    float clha = (sin(h) - sin(lat) * sin(moonDec)) / (cos(lat) * cos(moonDec));
    if (abs(clha) <= 1) clha = acos(clha);
    else
    {
        clha=0;
        /*Yet another thing to consider: the Sun is always in the south near 12:00 local time, but the Moon may be in the
         * south at any time of the day (or night). This means you must pay more attention that you're really iterating
         * towards the rise or set time you want, and not some rise/set time one day earlier, or later.

        Since the Moon rises and sets on the average 50 minutes later each day, there usually will be one day each month when
        the Moon never rises, and another day when it never sets. You must have this in mind when iterating your rise/set times,
        otherwise your program may easily get caught into an infinite loop when it tries to force e.g. a rise time between
        00:00 and 24:00 local time on a day when the Moon never rises.

        At high latitudes the Moon occasionally rises, or sets, twice on a single calendar day. This may happen above the
        "lunar arctic circle", which moves between 61.5 and 71.9 deg latitude during the 18-year period of the motion of
        the lunar nodes. You may want to pay attention to this.*/
    }
    moonRise = UT_Moon_in_south - clha * 12 / PI;
    moonSet = UT_Moon_in_south + clha * 12 / PI;
}

float Minutes(float bad_time){
 float good_time;
 int time_int;
 float help;
 float help_dec;
time_int=(int)bad_time;
help_dec=bad_time-time_int;
help=help_dec*60/100;
good_time=time_int+help;
  return good_time;
  }
void setup()
{
  Serial.begin(9600);
    int y=2018,m=4;
    for(int i=1; i<=30; i++)
    {
      Serial.println(String(i)+"."+String(m)+"."+String(y));

        sun(y, m, i, 12);
        moon(y, m, i, 12);
        if(sunRise>24){
          sunRise=sunRise-24;
          sunRise=Minutes(sunRise);
        Serial.println("Sun rise: "+String(sunRise)+"("+String(i+1)+"."+String(m)+"."+String(y)+")");
        }
        else if(sunRise<0){
          sunRise=24+sunRise;
          sunRise=Minutes(sunRise);
        Serial.println("Sun rise: "+String(sunRise)+"("+String(i-1)+"."+String(m)+"."+String(y)+")");
        }
        else {
          sunRise=Minutes(sunRise);
          Serial.println("Sun rise: "+String(sunRise));
          }
        if(sunSet>24){
          sunSet=sunSet-24;
          sunSet=Minutes(sunSet);
        Serial.println("Sun rise: "+String(sunSet)+"("+String(i+1)+"."+String(m)+"."+String(y)+")");
        }
        else if(sunSet<0){
          sunSet=24+sunSet;
            sunSet=Minutes(sunSet);
        Serial.println("Sun rise: "+String(sunSet)+"("+String(i-1)+"."+String(m)+"."+String(y)+")");
        }
        else {
            sunSet=Minutes(sunSet);
          Serial.println("Sun rise: "+String(sunSet));
        }



        float moonRise1=moonRise;
        moon(y,m,i,moonSet);
        float moonSet1=moonSet;
        moon(y,m,i,moonRise1);


      if(moonRise>24){
          moonRise=moonRise-24;
          moonRise=Minutes(moonRise);
        Serial.println("Moon rise: "+String(moonRise)+"("+String(i+1)+"."+String(m)+"."+String(y)+")");
        }
        else if(moonRise<0){
          moonRise=24+moonRise;
                    moonRise=Minutes(moonRise);
        Serial.println("Moon rise: "+String(moonRise)+"("+String(i-1)+"."+String(m)+"."+String(y)+")");
        }
        else {
          moonRise=Minutes(moonRise);
          Serial.println("Moon rise: "+String(moonRise));
        }
        if(moonSet1>24){

          moonSet1=moonSet1-24;
          moonSet1=Minutes(moonSet1);
        Serial.println("Moon set: "+String(moonSet1)+"("+String(i+1)+"."+String(m)+"."+String(y)+")");
        }
        else if(moonSet1<0){
          moonSet1=24+moonSet1;
            moonSet1=Minutes(moonSet1);
        Serial.println("Moon set: "+String(moonSet1)+"("+String(i-1)+"."+String(m)+"."+String(y)+")");
        }
        else {
          moonSet1=Minutes(moonSet1);
          Serial.println("Moon set: "+String(moonSet1));
          }

    }

}


void loop()
{
    // put your main code here, to run repeatedly:

}