/*  
 *  - - - - - - - - - -
 *   g a l _ m a f m s  
 *  - - - - - - - - - -
 *  
 *  This routine is part of the General Astrodynamics Library  
 *  
 *  Description: 
 *  
 *     This routine calculates the parameters for the Mars
 *     Fictitious Mean Sun and related parameters. 
 *  
 *  Status:
 *  
 *     support routine.  
 *  
 *  Given:
 *  
 *     tt1                 d        TT date part 1 ( see Note 1 )
 *     tt2                 d        TT date part 2 ( see Note 1 )  
 *  
 *  Returned:
 *  
 *     *m                  d        Mean anomaly ( radians )  
 *     *fms                d        Fictitious Maen Sun Angle ( radians )
 *     *pbs                d        Sum of angular perturbations in longitude ( radians )
 *     *ls                 d        Aerocentric solar longitude ( radians )
 *     *eot                d        Equation of Time ( radians ) 
 *  
 *  Notes:
 *
 *  1) The Julian Date is apportioned in any convenient way between  
 *     the arguments tt1 and tt2.  For example, JD=2450123.7 could  
 *     be expressed in any of these ways, among others:  
 *  
 *               tt1          tt2  
 *  
 *         2450123.7          0.0   (JD method)  
 *         2451545.0      -1421.3   (J2000 method)  
 *         2400000.5      50123.2   (MJD method)  
 *         2450123.5          0.2   (date & time method)  
 *
 *  Called:
 *
 *     gal_anp             Normalize angle to 0 <= a < pi
 *
 *  References:
 *
 *     A post-Pathfinder evaluation of areocentric solar coordinates with 
 *     improved timing recipes for Mars seasonal/diurnal climate studies 
 *     by Michael Allison, Megan McEwen, 
 *     Planetary and Space Science 48 (2000) 215-235 
 *
 *     Mars24 URL: http://www.giss.nasa.gov/tools/mars24/help/algorithm.html
 *     The referenced URL contains corrections to the referenced article.
 *
 *  This revision:
 *  
 *     2009 January 5  
 *  
 *  Copyright (C) 2009 Paul C. L. Willmott. See notes at end.  
 *  
 *-----------------------------------------------------------------------  
 */  
  
#include <math.h>
#include "gal_mafms.h"
#include "gal_const.h"
#include "gal_anp.h"
  
void 
gal_mafms  
 (  
    double tt1,
    double tt2,
    double *m,
    double *fms,
    double *pbs,
    double *ls,
    double *eot
 ) 
 
{

  double dt, x, p0, p1, p2, p3, p4, p5, p6, p7, vmm, a ;

  int i ;
  
/*
 * Primary short-term perturbations for the aerocentric solar longitude.
 * From Table 5 equivalent in URL referenced.
 */
 
  const double pert[7][3] = {
    { 0.0071,  2.2353,  49.409, } ,
    { 0.0057,  2.7543, 168.173, } ,
    { 0.0039,  1.1177, 191.837, } ,
    { 0.0037, 15.7866,  21.736, } ,
    { 0.0021,  2.1354,  15.704, } ,
    { 0.0020,  2.4694,  95.528, } ,
    { 0.0018, 32.8493,  49.095, } ,
  } ;

    
/*  
 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
 */  

/*
 * Calculate delta tt since J2000
 */
 
  dt = ( tt1 - GAL_J2000 ) + tt2 ;
   
/*
 * Determine Mars mean anomaly. (AM2000, eq. 16)
 */

  *m = gal_anp ( GAL_D2R * ( 19.3870 + 0.52402075 * dt ) ) ;

/*
 * Determine angle of Fiction Mean Sun. (AM2000, eq. 17)
 */
 
  *fms = gal_anp ( GAL_D2R * (270.3863 + 0.52403840 * dt ) ) ;

/*
 * Determine perturbers. (AM2000, eq. 18)
 */

  *pbs = 0.0 ;
   
  for ( i = 0; i < 7; i++ ) {
    a   = GAL_D2R * pert[i][0] ;
    x = GAL_D2R * ( 360.0 / 365.25 * dt / pert[i][1] + pert[i][2] ) ;
    *pbs += a * cos ( x ) ;
  } 

/*
 * Determine Equation of Center. (Bracketed term in AM2000, eqs. 19 and 20)
 * The equation of center is the true anomaly minus mean anomaly.
 */

  p0 = GAL_D2R * ( 10.691 + 3.0e-7 * dt ) ; 
  p1 = GAL_D2R * 0.623 ;
  p2 = GAL_D2R * 0.050 ;
  p3 = GAL_D2R * 0.005 ;
  p4 = GAL_D2R * 0.0005 ;
  
  vmm = p0 * sin ( (*m) ) + 
        p1 * sin ( 2.0 * (*m) ) + 
        p2 * sin ( 3.0 * (*m) ) + 
        p3 * sin ( 4.0 * (*m) ) + 
        p4 * sin ( 5.0 * (*m) ) + 
        (*pbs) ;

/*
 * Determine areocentric solar longitude. (AM2000, eq. 19)
 */
 
  *ls = gal_anp ( *fms + vmm ) ;
 
/*
 * Determine Equation of Time. (AM2000, eq. 20)
 */

  p5 = GAL_D2R * 2.861 ;
  p6 = GAL_D2R * 0.071 ;
  p7 = GAL_D2R * 0.002 ;
  
  *eot = p5 * sin ( 2.0 * (*ls) ) - 
         p6 * sin ( 4.0 * (*ls) ) + 
         p7 * sin ( 6.0 * (*ls) ) - 
         vmm ;

/*  
 * Finished.  
 */  

}

/*  
 *  gal - General Astrodynamics Library  
 *  Copyright (C) 2009 Paul C. L. Willmott  
 *  
 *  This program is free software; you can redistribute it and/or modify  
 *  it under the terms of the GNU General Public License as published by  
 *  the Free Software Foundation; either version 2 of the License, or  
 *  (at your option) any later version.  
 *  
 *  This program is distributed in the hope that it will be useful,  
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of  
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the  
 *  GNU General Public License for more details.  
 *  
 *  You should have received a copy of the GNU General Public License along  
 *  with this program; if not, write to the Free Software Foundation, Inc.,  
 *  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.  
 *  
 *  Contact:  
 *  
 *  Paul Willmott  
 *  vp9mu@amsat.org  
 */