/* pb0r.c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C-language version adapted from IDL code: $SSW/gen/idl/solar/pb0r.pro * Andrew L. Stanger HAO/NCAR 25 March 1999 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Project : SOHO - CDS * * Name : PB0R() * * Purpose : To calculate the solar P, B0 angles and the semi-diameter. * * Explanation : Uses semi-rigorous formulae to calculate the solar P (position * angle of pole) and B0 (latitude of point at disk centre) angles * and also the semi-diameter of the solar disk at the date/time * requested. * * Use : IDL> ang = pb0r(date_time,/soho) * * Inputs : date_time - the date/time specified in any CDS format * * Opt. Inputs : None * * Outputs : Function returns a 3-element array with * ang(0) = P (degrees) * ang(1) = B0 (degrees) * ang(2) = R semi-diameter (arcmin or * arcsec if * keyword set) * * Opt. Outputs: None * * Keywords : SOHO - if present the semi-diameter returned is as viewed * from the SOHO spacecraft * ARCSEC - returns semi-diameter in arcsec rather than arcmins * ERROR - Output keyword containing error message; * a null string is returned if no error occurs * * Calls : ANYTIM2UTC * SUN_POS * * Common : pb0r_common (internal common block) * * Restrictions: None * * Side effects: None * * Category : Util, coords * * Prev. Hist. : Based on Fortran programs by Hohenkerk and Emerson (RGO) * * Written : CDS/IDL version, C D Pike, RAL, 16-May-94 * * Modified : Update semi-diameter calculation, CDP, 20-May-94 * Version 3, William Thompson, GSFC, 14 November 1994 * Modified .DAY to .MJD * Version 4, CDP, 10-Jan-96 * Add SOHO/ARCSEC keywords and make NOW the default. * Version 5, Liyun Wang, GSFC/ARC, March 12, 1996 * Modified such that point of view can be changed to * SOHO if the env variable SC_VIEW is set to 1 * (via the call to USE_SOHO_VIEW) * Added ERROR keyword * Version 6, Liyun Wang, GSFC/ARC, March 14, 1996 * Replaced call to GET_ORBIT with the IDL call_function * Version 7, Liyun Wang, GSFC/ARC, March 21, 1996 * Modified such that if no orbit file is found, earth * view is used * Version 8, Liyun Wang, GSFC/ARC, April 10, 1996 * Set SC_VIEW to 0 if no orbit files are found * Version 9, February 6, 1997, Liyun Wang, NASA/GSFC * Changed call to ANYTIM2JD instead of CDS2JD * * Version : Version 9, February 6, 1997 */ #include #include int pb0r (iyear, imonth, iday, ihour, iminute, isecond, pangle, bangle, rsun) int iyear, imonth, iday, ihour, iminute, isecond ; double *pangle, *bangle, *rsun ; { extern double moduloc () ; extern double julday () ; static double pi = 3.14159265 ; double raddeg ; double degrad ; double jd ; double de ; double longmed ; double ra, dec ; double appl ; double oblt ; double lambda ; double node ; double arg ; double t ; double r ; double d ; double mv, me, mm, mj ; double p1, p2, p3, p4, p5, p6, p7, p8, p9 ; double p10, p11, p12, p13, p14, p15, p16 ; degrad = 180.0 / pi ; raddeg = pi / 180.0 ; /*--------------------------------------------------------------------------- * number of Julian days since 2415020.0 *--------------------------------------------------------------------------- * jd = anytim2jd (utc) * jd.int = jd.int - 2415020L */ jd = julday (iyear, imonth, iday, ihour, iminute, isecond) ; jd -= 2415020.0L ; /*--------------------------------------------------------------------------- * ignoring difference between UT and ET then ... *--------------------------------------------------------------------------- * de = DOUBLE (jd.int + jd.frac) */ de = (double) jd ; /*--------------------------------------------------------------------------- * get the longitude of the sun etc. *--------------------------------------------------------------------------- * sun_pos, de, longmed, ra, dec, appl, oblt */ sun_pos (de, &longmed, &ra, &dec, &appl, &oblt) ; /*--------------------------------------------------------------------------- * form aberrated longitude *--------------------------------------------------------------------------- */ lambda = longmed - (20.5L / 3600.0L) ; /*--------------------------------------------------------------------------- * form longitude of ascending node of sun's equator on ecliptic *--------------------------------------------------------------------------- */ node = 73.666666L + (50.25L / 3600.0L) * ( (de / 365.25L) + 50.0L ) ; arg = lambda - node ; /*--------------------------------------------------------------------------- * calculate P, the p1osition angle of the pole. *--------------------------------------------------------------------------- */ *pangle = (atan (-tan (oblt * raddeg) * cos (appl * raddeg)) + atan ( -0.12722L * cos (arg * raddeg))) * degrad ; p1 = -tan (oblt * raddeg) ; p2 = cos (appl * raddeg) ; p3 = p1 * p2 ; p4 = atan (p3) ; p5 = cos (arg * raddeg) ; p6 = atan ( -0.12722L * p5) ; p7 = p4 + p6 ; p8 = p7 * degrad ; /*--------------------------------------------------------------------------- * ... and B0 the tilt of the axis. *--------------------------------------------------------------------------- */ *bangle = asin ( 0.12620L * sin (arg * raddeg) ) * degrad ; /*--------------------------------------------------------------------------- * ... and the semi-diameter * * * Form the mean anomalies of Venus(MV),Earth(ME),Mars(MM),Jup1iter(MJ) * and the mean elongation of the Moon from the Sun(D). * *--------------------------------------------------------------------------- */ t = de / 36525.0L ; /* * mv = 212.6L + moduloc ( (58517.80L * t), 360.0L ) ; */ p1 = 212.6L ; p2 = 58517.80L ; p3 = 360.0L ; mv = p1 + moduloc ((p2 * t), p3) ; /* * me = 358.476L + moduloc ( (35999.0498L * t), 360.0L ) ; */ p1 = 358.476L ; p2 = 35999.0498L ; p3 = 360.0L ; me = p1 + moduloc ((p2 * t), p3) ; /* * mm = 319.5L + moduloc ( (19139.86L * t), 360.0L ) ; */ p1 = 319.5L ; p2 = 19139.86L ; p3 = 360.0L ; mm = p1 + moduloc ((p2 * t), p3) ; /* * mj = 225.3L + moduloc ( ( 3034.69L * t), 360.0L ) ; */ p1 = 225.3L ; p2 = 3034.69L ; p3 = 360.0L ; mj = p1 + moduloc ((p2 * t), p3) ; /* * d = 350.7L + moduloc ( (445267.11L * t), 360.0L ) ; */ p1 = 350.7L ; p2 = 445267.11L ; p3 = 360.0L ; d = p1 + moduloc ((p2 * t), p3) ; /*--------------------------------------------------------------------------- * Form the geocentric distance (r) and solar radius (rsun). *--------------------------------------------------------------------------- */ /* * r = 1.000141L - (0.016748L - 0.0000418L * t) * cos (me * raddeg) * - 0.000140L * cos (2.0L * me * raddeg) * + 0.000016L * cos ((58.3L + 2.0L * mv - 2.0L * me) * raddeg) * + 0.000005L * cos ((209.1L + mv - me)*raddeg) * + 0.000005L * cos ((253.8L - 2.0L * mm + 2.0L * me) * raddeg) * + 0.000016L * cos (( 89.5L - mj + me) * raddeg) * + 0.000009L * cos ((357.1L - 2.0L * mj + 2.0L * me) * raddeg) * + 0.000031L * cos (d * raddeg) ; */ p1 = 1.000141L ; p2 = 0.016748L ; p3 = 0.0000418L ; p4 = 0.000140L ; p5 = 2.0L ; p6 = 0.000016L ; p7 = 58.3L ; p8 = 2.0L ; p9 = 0.000005L ; p10 = 209.1L ; p11 = 253.8L ; p12 = 89.5L ; p13 = 0.000009L ; p14 = 357.1L ; p15 = 0.000031L ; r = p1 - (p2 - p3 * t) * cos (me * raddeg) - p4 * cos (p5 * me * raddeg) + p6 * cos ((p7 + p8 * mv - p8 * me) * raddeg) + p9 * cos ((p10 + mv - me) * raddeg) + p9 * cos ((p11 - p8 * mm + p8 * me) * raddeg) + p6 * cos ((p12 - mj + me) * raddeg) + p13 * cos ((p14 - p8 * mj + p8 * me) * raddeg) + p15 * cos (d * raddeg) ; /* *rsun = (0.2665685L / r) * 3600.0L ; */ /* arcseconds. */ p1 = 0.2665685L ; p2 = 3600.0L ; *rsun = (p1 / r) * p2 ; return (0) ; }