/*
 * Copyright (c) 1993, 1999 Alexandre Wennmacher (wennmach@geo.Uni-Koeln.DE)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 * This product includes software developed by Alexandre Wennmacher.
 * 4. The name of Alexandre Wennmacher may not be used to endorse or promote
 *    products derived from this software without specific prior written
 *    permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY ALEXANDRE WENNMACHER AND CONTRIBUTORS ``AS IS''
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL ALEXANDRE WENNMACHER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include <math.h>
#include "const.h"

double
EqunofTime(double JulianDate, double *SunRa, double *SunDec)
{
    double n, L, g, lambda, epsilon;

    n = JulianDate - 2451545.0;
    L = 280.460 + 0.9856474*n;
    g = 357.528 + 0.9856003*n;

    while (L < 0.0) L += 360.0;
    while (L >= 360.0) L -= 360.0;
    while (g < 0.0) g += 360.0;
    while (g > 360.0) g -= 360.0;

    lambda = L + 1.915*sin(DEGTORAD*g) + 0.020*sin(2.0*DEGTORAD*g);
    epsilon = 23.439 - 0.0000004*n;
  
    *SunRa  = atan(cos(DEGTORAD*epsilon)*tan(DEGTORAD*lambda));
    *SunDec = asin(sin(DEGTORAD*epsilon)*sin(DEGTORAD*lambda));
    if (*SunRa < 0.0) *SunRa += TWO_PI;
  
    lambda *= DEGTORAD;
    while (sin(lambda)*sin(*SunRa) < 0.0 || cos(lambda)*cos(*SunRa) < 0.0) {
        *SunRa += PI_2;
        if (*SunRa >= TWO_PI) *SunRa -= TWO_PI;
    }
  
    return(DEGTORAD*L - *SunRa);
}