;+ ; NAME: ; NewcombSun ; PURPOSE: ; Calculates the true ecliptic coordinates (longitude, latitude and ; distance) of the Sun according to Newcomb's theory of solar motion. ; See O. Montenbruck, Practical Ephemeris Calculations, par. 4.1.2, p. 66. ; CATEGORY: ; smei/gen/idl/ephem ; CALLING SEQUENCE: FUNCTION NewcombSun, UT, $ Precision = Precision , $ Longitude = Longitude , $ Latitude = Latitude , $ Distance = Distance , $ diskradius = DiskRadius , $ parallax = Parallax , $ degrees = Degrees , $ heliocentric=Heliocentric ; INPUTS: ; UT array; type: standard time structure ; times (UT) ; OPTIONAL INPUTS: ; /Precision if set, perturbations by several planets are ; taken into account ; /Longitude if set, return ecl. longitude ; /Latitude if set, return ecl. latitude ; /Distance if set, return Sun-Earth distance ; (in neither of the tree keywords LONGITUDE,LATITUDE, ; and DISTANCE is set, all three quantities are returned) ; /Degrees if set, return angles in degrees (default: radians) ; /Heliocentric if set, the longitude and latitude of Earth are returned ; (180 degrees is added to the longitude, the latitude ; is multiplied by -1, the distance is the same) ; OUTPUTS: ; Result array; type: double ; cliptic coordinates of the Sun ; ecliptic longitude and latitude (degrees or radians) ; and Sun-Earth distance in AU ; If only one coordinate (longitude, latitude or distance) is ; requested then the structure of the output array the same as for T. ; If more than one coordinate is requested the output array has an ; additional leading dimension with 2 or 3 elements. The leading dimension ; contains longitude, latitude and distance, in that order. ; OPTIONAL OUTPUTS: ; Diskradius radius of Sun in arcsec ; Parallax horizontal parallax in arcsec ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, TimeGet, ToRadians, SyncDims ; SEE ALSO: ; jpl_eph, jpl_test ; RESTRICTIONS: ; The calculation of Diskradius and Parallax needs the Sun-Earth ; distance. Hence the keyword /Distance must be set or all three ; keywords /Distance, /Longitude and /Latitude should be omitted. ; PROCEDURE: ; DLP long period perturbation of solar mean longitude and mean ; anomaly ; G mean anomaly of the Sun ; LO = DLO+SLO/3600 = mean longitude of Sun ; DL difference between true and mean solar longitude according to ; two body problem ; The precision calculation includes perturbations by planets, oscillatory ; terms with amplitudes less than 8". ; ; The alternative way to calculate the geocentric location of the ; Sun or the heliocentric location of the Earth is to use jpl_eph ; (or equivalently, big_eph). In most of our IDL code NewcombSun ; calls, have been replaced by big_eph calls. ; ; The following two calls should give essentially the same result: ; p = NewcombSun(tt,/precision,/degrees) ; p = big_eph(tt,body='sun',center='earth',/degrees,/to_sphere,/to_ecliptic,/precess) ; ; A check with the Astronomical Almanac shows: ; (All times are at midnight on the day given) ; ; Normal Precision Almanac ; 1971 Jan 1 Lng 279.9114 279.91575 279 54' 57.0" (279.91583) ; Lat 0. 0.38" 0.42" ; Dist 0.9832947 0.9833006 0.9832998 ; 1974 Jan 1 Lng 280.1882 280.1874 280 11' 14.94" (280.1875) ; Lat 0. ; Dist ; 1979 Jan 1 Lng 279.9702 279.9708 279 58' 14.90" (279.9708) ; Lat 0. ; Dist ; 1989 Jan 1 Lng 280.5535 280.5525 280 33' 10.44" (280.5529) ; Lat 0. ; Dist ; MODIFICATION HISTORY: ; 1989, Paul Hick (MPAE,UCSD) ; 10/24/91, Tom Davidson : accuracy test ; AUG-1993, Paul Hick, extension of SunEclLong ; FEB-1998, Paul Hick (UCSD/CASS, pphick@ucsd.edu) ; modified to accept multi-dimensional input arrays ; added keyword /heliocentric ; replaced keyword /radians by keyword /degrees ;- InitVar, Precision , /key InitVar, Longitude , /key InitVar, Latitude , /key InitVar, Distance , /key InitVar, Heliocentric , /key IF Longitude+Latitude+Distance EQ 0 THEN BEGIN Longitude = 1 Latitude = 1 Distance = 1 ENDIF rpd = ToRadians(/degrees) ; pi/180 in double precision ;T = TimeGet(UT, /jepoch, /scalar) ; Julian epoch ;T = (T-1900d0)/100d0 ; Time since 0.5 Jan 1900 (Julian centuries) T = 1.0d0+TimeGet(UT, /djepoch, /scalar)/100.0d0; Time since 0.5 Jan 1900 (Julian centuries) DLP = (1.882d0-0.016d0*T)*sin( ( 57.24d0+150.27d0*T)*rpd ) $ + 6.400d0*sin( (231.19d0+ 20.20d0*T)*rpd ) $ + 0.266d0*sin( ( 31.80d0+119.00d0*T)*rpd ) G = 358.4758333d0 + (35999.0d0*T mod 360.0d0) + (T*(179.10d0-0.54d0*T)+DLP)/3600.0d0 IF Precision THEN BEGIN G2 = 212.45d0+58517.493d0*T G4 = 319.58d0+19139.977d0*T G5 = 225.28d0+ 3034.583d0*T $ +1300.0d0*sin( (133.775d0+39.804d0*T)*rpd )/3600.0d0 G6 = 175.60d0+ 1221.794d0*T D = 350.737486d0+445267.114217d0*T A = 296.104608d0+477198.849108d0*T U = 11.250889d0+483202.025150d0*T ENDIF IF Longitude THEN BEGIN ; LO = DLO+SLO DLO = 279.6966778d0 + (36000.0d0*T mod 360.0d0) SLO = T*(2768.13d0+1.089d0*T) + .202d0*sin( (315.6d0+893.3d0*T)*rpd ) + DLP DL = (6910.057d0-17.240d0*T)*sin( G*rpd) $ + ( 72.338d0- 0.361d0*T)*sin(2.0d0*G*rpd) $ + 1.054d0*sin(3.0d0*G*rpd) IF Precision THEN BEGIN DL = DL $ ; DL2 (Venus) + 4.838d0*cos( (299.102d0+ G2- G)*rpd ) $ + 0.116d0*cos( (148.900d0+2.0d0*G2- G)*rpd ) $ + 5.526d0*cos( (148.313d0+2.0d0*G2-2.0d0*G)*rpd ) $ + 2.497d0*cos( (315.943d0+2.0d0*G2-3.0d0*G)*rpd ) $ + 0.666d0*cos( (177.710d0+3.0d0*G2-3.0d0*G)*rpd ) $ + 1.559d0*cos( (345.253d0+3.0d0*G2-4.0d0*G)*rpd ) $ + 1.024d0*cos( (318.150d0+3.0d0*G2-5.0d0*G)*rpd ) $ + 0.210d0*cos( (206.200d0+4.0d0*G2-4.0d0*G)*rpd ) $ + 0.144d0*cos( (195.400d0+4.0d0*G2-5.0d0*G)*rpd ) $ + 0.152d0*cos( (343.800d0+4.0d0*G2-6.0d0*G)*rpd ) $ + 0.123d0*cos( (195.300d0+5.0d0*G2-7.0d0*G)*rpd ) $ + 0.154d0*cos( (359.600d0+5.0d0*G2-8.0d0*G)*rpd ) $ ; DL4 (Mars) + 0.273d0*cos( (217.700d0- G4+ G)*rpd ) $ + 2.043d0*cos( (343.888d0-2.0d0*G4+2.0d0*G)*rpd ) $ + 1.770d0*cos( (200.402d0-2.0d0*G4+ G)*rpd ) $ + 0.129d0*cos( (294.200d0-3.0d0*G4+3.0d0*G)*rpd ) $ + 0.425d0*cos( (338.880d0-3.0d0*G4+2.0d0*G)*rpd ) $ + 0.500d0*cos( (105.180d0-4.0d0*G4+3.0d0*G)*rpd ) $ + 0.585d0*cos( (334.060d0-4.0d0*G4+2.0d0*G)*rpd ) $ + 0.204d0*cos( (100.800d0-5.0d0*G4+3.0d0*G)*rpd ) $ + 0.154d0*cos( (227.400d0-6.0d0*G4+4.0d0*G)*rpd ) $ + 0.101d0*cos( ( 96.300d0-6.0d0*G4+3.0d0*G)*rpd ) $ + 0.106d0*cos( (222.700d0-7.0d0*G4+4.0d0*G)*rpd ) $ ; DL5 (Jupiter) + 0.163d0*cos( (198.600d0- G5+2.0d0*G)*rpd ) $ + 7.208d0*cos( (179.532d0- G5+ G)*rpd ) $ + 2.600d0*cos( (263.217d0- G5 )*rpd ) $ + 2.731d0*cos( ( 87.145d0-2.0d0*G5+2.0d0*G)*rpd ) $ + 1.610d0*cos( (109.493d0-2.0d0*G5+ G)*rpd ) $ + 0.164d0*cos( (170.500d0-3.0d0*G5+3.0d0*G)*rpd ) $ + 0.556d0*cos( ( 82.650d0-3.0d0*G5+2.0d0*G)*rpd ) $ + 0.210d0*cos( ( 98.500d0-3.0d0*G5+ G)*rpd ) $ ; DL6 (Saturn) + 0.419d0*cos( (100.580d0- G6+ G)*rpd ) $ + 0.320d0*cos( (269.460d0- G6 )*rpd ) $ + 0.108d0*cos( (290.600d0-2.0d0*G6+2.0d0*G)*rpd ) $ + 0.112d0*cos( (293.600d0-2.0d0*G6+ G)*rpd ) $ ; DLM (Moon) + 6.454d0*sin( D *rpd ) $ + 0.177d0*sin( (D+A)*rpd ) $ - 0.424d0*sin( (D-A)*rpd ) $ + 0.172d0*sin( (D-G)*rpd ) ENDIF RL = DLO+(SLO+DL)/3600.0d0 IF Heliocentric THEN LngSun = RL+180.0d0 ELSE LngSun = RL LngSun = ( (LngSun mod 360.0d0)+360.0d0 ) mod 360.0d0 LngSun = LngSun*(rpd/ToRadians(degrees=Degrees)) ENDIF IF Distance THEN BEGIN R0 = ( 0.00003057d0-0.00000015d0*T) $ + (-0.00727412d0+0.00001814d0*T)*cos( G*rpd ) $ + (-0.00009138d0+0.00000046d0*T)*cos( 2.0d0*G*rpd ) $ + (-0.00000145d0) *cos( 3.0d0*G*rpd ) if Precision then begin R1 = $ ; DR2 (Venus) + 2359.0d0*cos( (209.080d0+ G2- G)*rpd ) $ + 160.0d0*cos( ( 58.400d0+2.0d0*G2- G)*rpd ) $ + 6842.0d0*cos( ( 58.318d0+2.0d0*G2-2.0d0*G)*rpd ) $ + 869.0d0*cos( (226.700d0+2.0d0*G2-3.0d0*G)*rpd ) $ + 1045.0d0*cos( ( 87.570d0+3.0d0*G2-3.0d0*G)*rpd ) $ + 1497.0d0*cos( (255.250d0+3.0d0*G2-4.0d0*G)*rpd ) $ + 194.0d0*cos( ( 49.500d0+3.0d0*G2-5.0d0*G)*rpd ) $ + 376.0d0*cos( (116.280d0+4.0d0*G2-4.0d0*G)*rpd ) $ + 196.0d0*cos( (105.200d0+4.0d0*G2-5.0d0*G)*rpd ) $ + 163.0d0*cos( (145.400d0+5.0d0*G2-5.0d0*G)*rpd ) $ + 141.0d0*cos( (105.400d0+5.0d0*G2-7.0d0*G)*rpd ) $ ; DR4 (Mars) + 150.0d0*cos( (127.700d0- G4+ G)*rpd ) $ + 2057.0d0*cos( (253.828d0-2.0d0*G4+2.0d0*G)*rpd ) $ + 151.0d0*cos( (295.000d0-2.0d0*G4+ G)*rpd ) $ + 168.0d0*cos( (203.500d0-3.0d0*G4+3.0d0*G)*rpd ) $ + 215.0d0*cos( (249.000d0-3.0d0*G4+2.0d0*G)*rpd ) $ + 478.0d0*cos( ( 15.170d0-4.0d0*G4+3.0d0*G)*rpd ) $ + 105.0d0*cos( ( 65.900d0-4.0d0*G4+2.0d0*G)*rpd ) $ + 107.0d0*cos( (324.600d0-5.0d0*G4+4.0d0*G)*rpd ) $ + 139.0d0*cos( (137.300d0-6.0d0*G4+4.0d0*G)*rpd ) $ ; DR5 (Jupiter) + 208.0d0*cos( (112.000d0- G5+2.0d0*G)*rpd ) $ + 7067.0d0*cos( ( 89.545d0- G5+ G)*rpd ) $ + 244.0d0*cos( (338.600d0- G5 )*rpd ) $ + 103.0d0*cos( (350.500d0-2.0d0*G5+3.0d0*G)*rpd ) $ + 4026.0d0*cos( (357.108d0-2.0d0*G5+2.0d0*G)*rpd ) $ + 1459.0d0*cos( ( 19.467d0-2.0d0*G5+ G)*rpd ) $ + 281.0d0*cos( ( 81.200d0-3.0d0*G5+3.0d0*G)*rpd ) $ + 803.0d0*cos( (352.560d0-3.0d0*G5+2.0d0*G)*rpd ) $ + 174.0d0*cos( ( 8.600d0-3.0d0*G5+ G)*rpd ) $ + 113.0d0*cos( (347.700d0-4.0d0*G5+2.0d0*G)*rpd ) $ ; DR6 (Saturn) + 429.0d0*cos( ( 10.600d0- G6+ G)*rpd ) $ + 162.0d0*cos( (200.600d0-2.0d0*G6+2.0d0*G)*rpd ) $ + 112.0d0*cos( (203.100d0-2.0d0*G6+ G)*rpd ) $ ; DRM (Moon) + 13360.0d0*cos( (D )*rpd ) $ + 370.0d0*cos( (D+A)*rpd ) $ - 1330.0d0*cos( (D-A)*rpd ) $ - 140.0d0*cos( (D+G)*rpd ) $ + 360.0d0*cos( (D-G)*rpd ) R0 = R0+R1*1.0d-9 ENDIF DisSun = 10^R0 DiskRadius = asin(695990.0d0/(DisSun*149597870L))/rpd*3600 Parallax = 8.794148d0/DisSun ENDIF IF Latitude THEN BEGIN LatSun = 0.0d0*T IF Precision THEN LatSun = $ -0.210d0*cos( (151.800d0+3.0d0*G2-4.0d0*G)*rpd ) $ -0.166d0*cos( (265.500d0-2.0d0*G5+ G)*rpd ) $ +0.576d0*sin( U*rpd ) ; Latitude in arcsec IF Heliocentric THEN LatSun = -LatSun LatSun = LatSun/3600*(rpd/ToRadians(degrees=Degrees)) ;LatSun = LatSun*(rpd/ToRadians(degrees=Degrees)) ENDIF n = Longitude+Latitude+Distance SyncDims, sz, replicateinfo=[size(indgen(n)),size(T)] Out = reform( replicate(0.0d0,sz[sz[0]+2]), n, sz[sz[0]+2]/n ) i = 0 IF Longitude THEN BEGIN Out[i,*] = LngSun & i = i+1 & ENDIF IF Latitude THEN BEGIN Out[i,*] = LatSun & i = i+1 & ENDIF IF Distance THEN Out[i,*] = DisSun Out = reform( Out, sz[1:sz[0]] ) IF n EQ 1 THEN SyncDims, Out, size=size(T) RETURN, Out & END