/* Copyright (GPL) 2004   Mike Chirico mchirico@comcast.net
   Updated: Sun Nov 28 15:15:05 EST 2004

   Program adapted by Mike Chirico mchirico@comcast.net

   Reference:
    http://prdownloads.sourceforge.net/souptonuts/working_with_time.tar.gz?download
    http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html


  Compile:

     gcc -o sunrise -Wall -W -O2 -s -pipe -lm sunrise.c

  or for debug output

     gcc -o sunrise -DDEBUG=1 -Wall -W -O2 -s -pipe -lm sunrise.c




  This can also go in a batch job to calculate the next
  20 days as follows:

    #!/bin/bash
    lat=39.95
    long=75.15
    for (( i=0; i <= 20; i++))
    do
     ./sunrise    `date -d "+$i day" "+%Y %m %d"` $lat $long
    done


   
*/

/* gcc -DDEBUG=1 .. */
#ifndef DEBUG
#define DEBUG 0
#endif


#include <sys>
#include <time>
#include <stdlib>
#include <stdio>
#include <math>






double calcSunEqOfCenter(double t);


/* Convert degree angle to radians */

double  degToRad(double angleDeg)
{
  return (M_PI * angleDeg / 180.0);
}

double radToDeg(double angleRad)
{
  return (180.0 * angleRad / M_PI);
}


double calcMeanObliquityOfEcliptic(double t)
{
  double seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
  double e0 = 23.0 + (26.0 + (seconds/60.0))/60.0;

  return e0;              // in degrees
}


double calcGeomMeanLongSun(double t)
{


  double L = 280.46646 + t * (36000.76983 + 0.0003032 * t);
  while( (int) L >  360 )
    {
      L -= 360.0;

    }
  while(  L <  0)
    {
      L += 360.0;

    }


  return L;              // in degrees
}






double calcObliquityCorrection(double t)
{
  double e0 = calcMeanObliquityOfEcliptic(t);


  double omega = 125.04 - 1934.136 * t;
  double e = e0 + 0.00256 * cos(degToRad(omega));
  return e;               // in degrees
}

double calcEccentricityEarthOrbit(double t)
{
  double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  return e;               // unitless
}

double calcGeomMeanAnomalySun(double t)
{
  double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
  return M;               // in degrees
}


double calcEquationOfTime(double t)
{


  double epsilon = calcObliquityCorrection(t);               
  double  l0 = calcGeomMeanLongSun(t);
  double e = calcEccentricityEarthOrbit(t);
  double m = calcGeomMeanAnomalySun(t);
  double y = tan(degToRad(epsilon)/2.0);
  y *= y;
  double sin2l0 = sin(2.0 * degToRad(l0));
  double sinm   = sin(degToRad(m));
  double cos2l0 = cos(2.0 * degToRad(l0));
  double sin4l0 = sin(4.0 * degToRad(l0));
  double sin2m  = sin(2.0 * degToRad(m));
  double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
            - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;


  return radToDeg(Etime)*4.0;   // in minutes of time


   }


double calcTimeJulianCent(double jd)
{

  double T = ( jd - 2451545.0)/36525.0;
  return T;
}









double calcSunTrueLong(double t)
{
  double l0 = calcGeomMeanLongSun(t);
  double c = calcSunEqOfCenter(t);

  double O = l0 + c;
  return O;               // in degrees
}



double calcSunApparentLong(double t)
{
  double o = calcSunTrueLong(t);

  double  omega = 125.04 - 1934.136 * t;
  double  lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega));
  return lambda;          // in degrees
}




double calcSunDeclination(double t)
{
  double e = calcObliquityCorrection(t);
  double lambda = calcSunApparentLong(t);

  double sint = sin(degToRad(e)) * sin(degToRad(lambda));
  double theta = radToDeg(asin(sint));
  return theta;           // in degrees
}


double calcHourAngleSunrise(double lat, double solarDec)
{
  double latRad = degToRad(lat);
  double sdRad  = degToRad(solarDec);



  double HA = (acos(cos(degToRad(90.833))/(cos(latRad)*cos(sdRad))-tan(latRad) * tan(sdRad)));

  return HA;              // in radians
}

double calcHourAngleSunset(double lat, double solarDec)
{
  double latRad = degToRad(lat);
  double sdRad  = degToRad(solarDec);


  double HA = (acos(cos(degToRad(90.833))/(cos(latRad)*cos(sdRad))-tan(latRad) * tan(sdRad)));

  return -HA;              // in radians
}


double calcJD(int year,int month,int day)
   {
      if (month <2>tm_hour;

  if(DEBUG)
    printf("delta=%d dst=%d\n",delta,dst);
 
  tseconds= seconds;
  if(DEBUG)
   printf("Number of seconds %ld\n",seconds);
  seconds= seconds + calcSunriseUTC( JD,  latitude,  longitude)*60;
  seconds= seconds - delta*3600;




  strftime(buffer,30,"%m-%d-%Y  %T",localtime(&seconds));
  printf("Sunrise  %s   ",buffer);



  seconds=tseconds;
  seconds+=calcSunsetUTC( JD,  latitude,  longitude)*60;
  seconds= seconds - delta*3600;


  strftime(buffer,30,"%m-%d-%Y  %T",localtime(&seconds));
  printf("Sunset %s\n",buffer);




  return 0;

}

#! /bin/bash
# /root/timeup 
# Author: fvillant
# Date: 2007.10.10.
# Recuperation du lever et du coucher du Soleil a Porto-Vecchio pour la date du lendemain

/root/sunrise `date -d "+1 day" "+%Y %m %d"` 41.41 -9.22

#!/usr/bin/perl -w
#
# 28 Janvier 2010 Creation Flavien Villant
# Edition crontab pour modification sunset
# sunrise selon le jour et l'heure la la
#
#
#
#
system("/root/timeup > sun.txt");
open(FILE,"sun.txt");
@arr = <FILE>;
for $line (@arr) {
chomp $line;
# Recuperation des horaires
$sunriseH=substr($line,21,2);
$sunriseM=substr($line,24,2);
$sunsetH=substr($line,51,2);
$sunsetM=substr($line,54,2);



$sunsetrise = "# Mode Jour/Nuit Zoneminder ".$Line." NE PAS DECOMMENTER";
$jour = $sunriseM." ".$sunriseH." * * * root /usr/bin/zmpkg.pl Day";
$nuit = $sunsetM." ".$sunsetH." * * * root /usr/bin/zmpkg.pl Night";

}

close(FILE);




open(FILE,"/etc/crontab");
my @crn = <FILE>;
close(FILE);

$crn[23] = $sunsetrise;
$crn[24] = $jour;
$crn[25] = $nuit;

open(FILE, ">/etc/crontab") || die "Erreur E/S:$!\n";
foreach my $ligne (@crn) {
    chomp $ligne;
    print FILE "$ligne\n";
}


system("chmod 754 /etc/crontab");
system("chmod +x /etc/crontab");

exit;