/* sun_pos.c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C-language version adapted from IDL code: $SSW/gen/idl/solar/sun_pos.pro * Andrew L. Stanger HAO/NCAR 25 March 1999 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Project : SOHO - CDS * * Name : SUN_POS * * Purpose : Calculate solar ephemeris parameters. * * Explanation : Allows for planetary and lunar perturbations in the calculation * of solar longitude at date and various other solar positional * parameters. * * Use : sun_pos (date, longitude, ra, dec, app_long, obliq) ; * * Inputs : date - fractional number of days since JD 2415020.0 * * Opt. Inputs : None * * Outputs : longitude - Longitude of sun for mean equinox of date (degs) * ra - Apparent RA for true equinox of date (degs) * dec - Apparent declination for true equinox of date (degs) * app_long - Apparent longitude (degs) * obliq - True obliquity (degs) * * Opt. Outputs: All above * * Keywords : None * * Calls : None * * Common : None * * Restrictions: None * * Side effects: None * * Category : Util, coords * * Prev. Hist. : From Fortran routine by B Emerson (RGO). * * Written : CDS/IDL version by C D Pike, RAL, 17-May-94 * * Modified : * * Version : Version 1, 17-May-94 * C-language version: Andrew L. Stanger, HAO/NCAR, 26 March 1999. */ /* pro sun_pos, dd, longmed, ra, dec, l, oblt */ #include int sun_pos (jd, longmed, ra, dec, app_long, obliq) double jd ; double *longmed ; double *ra ; double *dec ; double *app_long ; double *obliq ; { extern double moduloc () ; double degrad ; double raddeg ; static double pi = 3.14159265 ; double d ; double ellcor ; double jupcorr ; double l ; double longterm ; double me ; double mj ; double mm ; double mv ; double marscorr ; double mooncorr ; double omega ; double vencorr ; double t ; double p1, p2, p3, p4, p5, p6, p7, p8, p9 ; double p10, p11, p12, p13, p14 ; double mp ; degrad = 180.0 / pi ; raddeg = pi / 180.0 ; /* * This routine is a truncated version of Newcomb's Sun and * is designed to give apparent angular coordinates (T.E.D) to a * precision of one second of time. */ /* Form time in Julian centuries from 1900.0 */ t = jd / 36525.0L ; /* Form sun's mean longitude. */ /* * l = (279.696678L + moduloc ((36000.768925L * t), 360.0L)) * 3600.0L ; */ p1 = 279.696678L ; p2 = 36000.768925L ; p3 = 360.0L ; p4 = 3600.0L ; l = (p1 + moduloc ((p2 * t), p3)) * p4 ; /* Allow for ellipticity of the orbit (equation of centre) * using the Earth's mean anomoly ME. */ /* * me = 358.475844L + moduloc ((35999.049750L * t), 360.0L) ; * ellcor = (6910.1L - 17.2L * t) * sin (me * raddeg) + * 72.3L * sin (2.0L * me * raddeg) ; */ p1 = 358.475844L ; p2 = 35999.049750L ; p3 = 360.0L ; me = p1 + moduloc ((p2 * t), p3) ; p1 = 6910.1L ; p2 = 17.2L ; p3 = 72.3L ; p4 = 2.0L ; ellcor = (p1 - p2 * t) * sin (me * raddeg) + p3 * sin (p4 * me * raddeg) ; l = l + ellcor ; /* Allow for the Venus perturbations using the mean anomaly of Venus MV. */ /* * mv = 212.603219L + moduloc ((58517.803875L * t), 360.0L) ; * vencorr = 4.8L * cos ((299.1017L + mv - me) * raddeg) + * 5.5L * cos ((148.3133L + 2.0L * mv - 2.0L * me ) * raddeg) + * 2.5L * cos ((315.9433L + 2.0L * mv - 3.0L * me ) * raddeg) + * 1.6L * cos ((345.2533L + 3.0L * mv - 4.0L * me ) * raddeg) + * 1.0L * cos ((318.15L + 3.0L * mv - 5.0L * me ) * raddeg) ; */ p1 = 212.603219L ; p2 = 58517.803875L ; p3 = 360.0L ; mv = p1 + moduloc ((p2 * t), p3) ; p1 = 4.8L ; p2 = 299.1017L ; p3 = 5.5L ; p4 = 148.3133L ; p5 = 2.0L ; p6 = 2.5L ; p7 = 315.9433L ; p8 = 3.0L ; p9 = 1.6L ; p10 = 345.2533L ; p11 = 4.0L ; p12 = 1.0L ; p13 = 318.15L ; p14 = 5.0L ; vencorr = p1 * cos ((p2 + mv - me) * raddeg) + p3 * cos ((p4 + p5 * mv - p5 * me) * raddeg) + p6 * cos ((p7 + p5 * mv - p8 * me) * raddeg) + p9 * cos ((p10 + p8 * mv - p11 * me) * raddeg) + p12 * cos ((p13 + p8 * mv - p14 * me) * raddeg) ; l = l + vencorr ; /* Allow for the Mars perturbations using the mean anomaly of Mars MM. */ /* * mm = 319.529425L + moduloc (( 19139.858500L * t), 360.0L ) ; */ p1 = 319.529425L ; p2 = 19139.858500L ; p3 = 360.0L ; mm = p1 + moduloc ((p2 * t), p3) ; /* * marscorr = 2.0L * cos ((343.8883L - 2.0L * mm + 2.0L * me) * raddeg ) + * 1.8L * cos ((200.4017L - 2.0L * mm + me) * raddeg) ; */ p1 = 2.0L ; p2 = 343.8883L ; p3 = 1.8L ; p4 = 200.4017L ; marscorr = p1 * cos ((p2 - p1 * mm + p1 * me) * raddeg) + p3 * cos ((p4 - p1 * mm + me) * raddeg) ; l = l + marscorr ; /* Allow for the Jupiter perturbations * using the mean anomaly of Jupiter mj. */ /* * mj = 225.328328L + moduloc (( 3034.6920239L * t), 360.0L ) ; */ p1 = 225.328328L ; p2 = 3034.6920239L ; p3 = 360.0L ; mj = p1 + moduloc ((p2 * t), p3) ; /* * jupcorr = 7.2L * cos (( 179.5317L - mj + me ) * raddeg) + * 2.6L * cos ((263.2167L - mj ) * raddeg) + * 2.7L * cos (( 87.1450L - 2.0L * mj + 2.0L * me ) * raddeg) + * 1.6L * cos ((109.4933L - 2.0L * mj + me ) * raddeg) ; */ p1 = 7.2L ; p2 = 179.5317L ; p3 = 2.6L ; p4 = 263.2167L ; p5 = 2.7L ; p6 = 87.1450L ; p7 = 2.0L ; p8 = 1.6L ; p9 = 109.4933L ; jupcorr = p1 * cos ((p2 - mj + me) * raddeg) + p3 * cos ((p4 - mj) * raddeg) + p5 * cos ((p6 - p7 * mj + p7 * me) * raddeg) + p8 * cos ((p9 - p7 * mj + me) * raddeg) ; l = l + jupcorr ; /* Allow for the Moons perturbations using the mean elongation of * the Moon from the Sun D. */ /* * d = 350.7376814L + moduloc (( 445267.11422L * t), 360.0L ) ; */ p1 = 350.7376814L ; p2 = 445267.11422L ; p3 = 360.0L ; d = p1 + moduloc ((p2 * t), p3) ; /* * mooncorr = 6.5L * sin (d * raddeg) ; */ p1 = 6.5L ; mooncorr = p1 * sin (d * raddeg) ; l = l + mooncorr ; /* Allow for long period terms. */ /* * longterm = + 6.4L * sin (( 231.19L + 20.20L * t ) * raddeg) ; */ p1 = 6.4L ; p2 = 231.19L ; p3 = 20.20L ; longterm = p1 * sin ((p2 + p3 * t) * raddeg) ; l = l + longterm ; /* * l = moduloc ((l + 2592000.0L), 1296000.0L) ; */ p1 = 2592000.0L ; p2 = 1296000.0L ; l = moduloc ((l + p1), p2) ; *longmed = l / 3600.0L ; /* Allow for Aberration. */ /* * l = l - 20.5L ; */ p1 = 20.5L ; l = l - p1 ; /* Allow for Nutation using the longitude of the Moons mean node OMEGA. */ /* * omega = 259.183275L - moduloc (( 1934.142008L * t ), 360.0L ) ; */ p1 = 259.183275L ; p2 = 1934.142008L ; p3 = 360.0L ; omega = p1 - moduloc ((p2 * t), p3) ; l = l - 17.2L * sin (omega * raddeg) ; /* Form the True Obliquity. */ /* * *obliq = 23.452294L - 0.0130125L * t + * (9.2L * cos (omega * raddeg)) / 3600.0L ; */ p1 = 23.452294L ; p2 = 0.0130125L ; p3 = 9.2L ; p4 = 3600.0L ; *obliq = p1 - p2 * t + (p3 * cos (omega * raddeg)) / p4 ; /* Form Right Ascension and Declination. */ /* * *app_long = l / 3600.0L ; */ p1 = 3600.0L ; *app_long = l / p1 ; /* * l = l / 3600.0L ; */ p1 = 3600.0L ; l = l / p1 ; *ra = atan2 ( sin (l * raddeg) * cos (*obliq * raddeg), cos (l * raddeg) ) * degrad ; if (*ra < 0.0L) *ra = *ra + 360.0L ; *dec = asin (sin (l * raddeg) * sin (*obliq * raddeg)) * degrad ; return (0) ; }