C+
C NAME:
C	SunNewcomb8
C PURPOSE:
C	Calculates the true ecliptic coordinates (longitude, latitude and
C	distance) of the Sun according to Newcomb's theory of solar motion.
C	See O. Montenbruck, Practical Ephemeris Calculations, par. 4.1.2, p. 66.
C CATEGORY:
C	Celestial mechanics
C CALLING SEQUENCE:
	subroutine SunNewcomb8(ID,iYr,Doy8,rLng,rLat,rDis)
C CALLS:
C	Julian
C INPUTS:
C	ID	integer		0 = return two-body solution
C				1 = take planetary perturbations into account
C	iYr	integer		year
C	Doy8	real*8		doy8 of year, including fraction for time of day
C OUTPUTS:
C	rLng	real*8	ecliptic longitude in degrees [0,360)
C	rLat	real*8	ecliptic latitude  in arcsec
C	rDis	real*8	Sun-Earth distance in AU
C PROCEDURE:
C	DLP	long period perturbation of solar mean longitude and mean
C		anomaly
C	G	mean anomaly of the Sun
C	LO	= DLO+SLO/3600 = mean longitude of Sun
C	DL	difference between true and mean solar longitude according to
C		two body problem
C	The precision calculation includes perturbations by planets, oscillatory
C	terms with amplitudes less than 8".
C ACCURACY:
C	A check with the Astronomical Almanac shows:
C	(All times are at midnight on the day given)
C
C	  			Normal		Precision	Almanac
C	1971	Jan 1	Lng	279.9114	279.91575	279 54' 57.0"  (279.91583)
C			Lat	  0.		0.38"		0.42"
C			Dist	0.9832947	0.9833006	0.9832998
C	1974	Jan 1	Lng	280.1882	280.1874	280 11' 14.94" (280.1875)
C			Lat	  0.
C			Dist
C	1979	Jan 1	Lng	279.9702	279.9708	279 58' 14.90" (279.9708)
C			Lat	  0.
C			Dist
C	1989	Jan 1	Lng	280.5535	280.5525	280 33' 10.44" (280.5529) 
C			Lat	  0.
C			Dist
C MODIFICATION HISTORY:
C	1989, Paul Hick (MPAE,UCSD/CASS; pphick@ucsd.edu)
C	10/24/91, Tom Davidson : accuracy test
C	AUG-1993, Paul Hick, extension of SunEclLong
C-
	    integer	ID
	    integer	iYr
	    real*8	Doy8
	    real*8	rLng
	    real*8	rLat
	    real*8	rDis

	real*8		JD
	real*8		T
	real*8		DLP
	real*8		G
	real*8		G2
	real*8		G4
	real*8		G5
	real*8		G6
	real*8		D
	real*8		A
	real*8		U
	real*8		DLO
	real*8		SLO
	real*8		DL
	real*8		R1

	call Julian8(0,iYr,Doy8,JD,T)
	T = (T-1900)/100	! Time since 0.5 Jan 1900 (Julian centuries)

	DLP =   (1.882d0-.016d0*T)*dsind( ( 57.24d0+150.27d0*T) )
     &		+ 6.40d0*dsind( (231.19d0+ 20.20d0*T) )
     &		+ .266d0*dsind( ( 31.80d0+119.00d0*T) )

	G = 358.4758333d0 + mod(35999.d0*T,360.d0) + (T*(179.10d0-.54d0*T)+DLP)/3600.d0

	!-------
	! LO = DLO+SLO

	DLO = 279.6966778d0 + mod(36000.d0*T,360.d0)
	SLO = T*(2768.13d0+1.089d0*T) + .202d0*dsind( 315.6d0+893.3d0*T ) + DLP

	DL =	  (6910.057d0-17.240d0*T)*dsind(     G)
     &		+ (  72.338d0- 0.361d0*T)*dsind(2.d0*G)
     &				+ 1.054d0*dsind(3.d0*G)

	rDis =( .00003057d0-.00000015d0*T)
     &		+ (-.00727412d0+.00001814d0*T)*dcosd(      G )
     &		+ (-.00009138d0+.00000046d0*T)*dcosd( 2.d0*G )
     &		+ (-.00000145d0)	      *dcosd( 3.d0*G )

	rLat = 0.d0

	if (ID .eq. 1) then

	    G2 = 212.45d0+58517.493d0*T

	    G4 = 319.58d0+19139.977d0*T
	    G5 = 225.28d0+ 3034.583d0*T	
     &		+ 1300.d0*dsind( 133.775d0+39.804d0*T )/3600.d0
	    G6 = 175.60d0+ 1221.794d0*T

	    D = 350.737486d0+445267.114217d0*T
	    A = 296.104608d0+477198.849108d0*T
	    U =  11.250889d0+483202.025150d0*T

	    !-------
	    ! DL2 (Venus)

	    DL = DL
     &		+4.838d0*dcosd( 299.102d0+     G2-     G )
     &		+0.116d0*dcosd( 148.900d0+2.d0*G2-     G )
     &		+5.526d0*dcosd( 148.313d0+2.d0*G2-2.d0*G )
     &		+2.497d0*dcosd( 315.943d0+2.d0*G2-3.d0*G )
     &		+0.666d0*dcosd( 177.710d0+3.d0*G2-3.d0*G )
     &		+1.559d0*dcosd( 345.253d0+3.d0*G2-4.d0*G )
     &		+1.024d0*dcosd( 318.150d0+3.d0*G2-5.d0*G )
     &		+0.210d0*dcosd( 206.200d0+4.d0*G2-4.d0*G )
     &		+0.144d0*dcosd( 195.400d0+4.d0*G2-5.d0*G )
     &		+0.152d0*dcosd( 343.800d0+4.d0*G2-6.d0*G )
     &		+0.123d0*dcosd( 195.300d0+5.d0*G2-7.d0*G )
     &		+0.154d0*dcosd( 359.600d0+5.d0*G2-8.d0*G )

	    !-------
	    ! DL4 (Mars)

	    DL = DL
     &		+0.273d0*dcosd( 217.700d0-     G4+     G )
     &		+2.043d0*dcosd( 343.888d0-2.d0*G4+2.d0*G )
     &		+1.770d0*dcosd( 200.402d0-2.d0*G4+     G )
     &		+0.129d0*dcosd( 294.200d0-3.d0*G4+3.d0*G )
     &		+0.425d0*dcosd( 338.880d0-3.d0*G4+2.d0*G )
     &		+0.500d0*dcosd( 105.180d0-4.d0*G4+3.d0*G )
     &		+0.585d0*dcosd( 334.060d0-4.d0*G4+2.d0*G )
     &		+0.204d0*dcosd( 100.800d0-5.d0*G4+3.d0*G )
     &		+0.154d0*dcosd( 227.400d0-6.d0*G4+4.d0*G )
     &		+0.101d0*dcosd(  96.300d0-6.d0*G4+3.d0*G )
     &		+0.106d0*dcosd( 222.700d0-7.d0*G4+4.d0*G )

	    !-------
	    ! DL5 (Jupiter)

	    DL = DL
     &		+0.163d0*dcosd( 198.600d0-     G5+2.d0*G )
     &		+7.208d0*dcosd( 179.532d0-     G5+     G )
     &		+2.600d0*dcosd( 263.217d0-     G5        )
     &		+2.731d0*dcosd(  87.145d0-2.d0*G5+2.d0*G )
     &		+1.610d0*dcosd( 109.493d0-2.d0*G5+     G )
     &		+0.164d0*dcosd( 170.500d0-3.d0*G5+3.d0*G )
     &		+0.556d0*dcosd(  82.650d0-3.d0*G5+2.d0*G )
     &		+0.210d0*dcosd(  98.500d0-3.d0*G5+     G )

	    !-------
	    ! DL6 (Saturn)

	    DL = DL
     &		+0.419d0*dcosd( 100.580d0-     G6+     G )
     &		+0.320d0*dcosd( 269.460d0-     G6        )
     &		+0.108d0*dcosd( 290.600d0-2.d0*G6+2.d0*G )
     &		+0.112d0*dcosd( 293.600d0-2.d0*G6+     G )

	    !-------
	    ! DLM (Moon)

	    DL = DL
     &		+6.454d0*dsind( D   )
     &		+0.177d0*dsind( D+A )
     &		-0.424d0*dsind( D-A )
     &		+0.172d0*dsind( D-G )

	    !-------
	    ! DR2 (Venus)

	    R1  =
     &		 2359.d0*dcosd( 209.080d0+     G2-     G )
     &		+ 160.d0*dcosd(  58.400d0+2.d0*G2-     G )
     &		+6842.d0*dcosd(  58.318d0+2.d0*G2-2.d0*G )
     &		+ 869.d0*dcosd( 226.700d0+2.d0*G2-3.d0*G )
     &		+1045.d0*dcosd(  87.570d0+3.d0*G2-3.d0*G )
     &		+1497.d0*dcosd( 255.250d0+3.d0*G2-4.d0*G )
     &		+ 194.d0*dcosd(  49.500d0+3.d0*G2-5.d0*G )
     &		+ 376.d0*dcosd( 116.280d0+4.d0*G2-4.d0*G )
     &		+ 196.d0*dcosd( 105.200d0+4.d0*G2-5.d0*G )
     &		+ 163.d0*dcosd( 145.400d0+5.d0*G2-5.d0*G )
     &		+ 141.d0*dcosd( 105.400d0+5.d0*G2-7.d0*G )

	    !-------
	    ! DR4 (Mars)

	    R1 = R1
     &		+ 150.d0*dcosd( 127.700d0-     G4+     G )
     &		+2057.d0*dcosd( 253.828d0-2.d0*G4+2.d0*G )
     &		+ 151.d0*dcosd( 295.000d0-2.d0*G4+     G )
     &		+ 168.d0*dcosd( 203.500d0-3.d0*G4+3.d0*G )
     &		+ 215.d0*dcosd( 249.000d0-3.d0*G4+2.d0*G )
     &		+ 478.d0*dcosd(  15.170d0-4.d0*G4+3.d0*G )
     &		+ 105.d0*dcosd(  65.900d0-4.d0*G4+2.d0*G )
     &		+ 107.d0*dcosd( 324.600d0-5.d0*G4+4.d0*G )
     &		+ 139.d0*dcosd( 137.300d0-6.d0*G4+4.d0*G )

	    !-------
	    ! DR5 (Jupiter)

	    R1 = R1
     &		+ 208.d0*dcosd( 112.000d0-     G5+2.d0*G )
     &		+7067.d0*dcosd(  89.545d0-     G5+     G )
     &		+ 244.d0*dcosd( 338.600d0-     G5        )
     &		+ 103.d0*dcosd( 350.500d0-2.d0*G5+3.d0*G )
     &		+4026.d0*dcosd( 357.108d0-2.d0*G5+2.d0*G )
     &		+1459.d0*dcosd(  19.467d0-2.d0*G5+     G )
     &		+ 281.d0*dcosd(  81.200d0-3.d0*G5+3.d0*G )
     &		+ 803.d0*dcosd( 352.560d0-3.d0*G5+2.d0*G )
     &		+ 174.d0*dcosd(   8.600d0-3.d0*G5+     G )
     &		+ 113.d0*dcosd( 347.700d0-4.d0*G5+2.d0*G )

	    !-------
	    ! DR6 (Saturn)

	    R1 = R1
     &		+ 429.d0*dcosd(  10.600d0-     G6+     G )
     &		+ 162.d0*dcosd( 200.600d0-2.d0*G6+2.d0*G )
     &		+ 112.d0*dcosd( 203.100d0-2.d0*G6+     G )

	    !-------
	    ! DRM (Moon)

	    R1 = R1
     &		+13360.d0*dcosd( D   )
     &		+  370.d0*dcosd( D+A )
     &		- 1330.d0*dcosd( D-A )
     &		-  140.d0*dcosd( D+G )
     &		+  360.d0*dcosd( D-G )

	    rDis = rDis+R1*1.d-9

	    rLat =
     &		-0.210d0*dcosd( 151.800d0+3.d0*G2-4.d0*G )
     &		-0.166d0*dcosd( 265.500d0-2.d0*G5+     G )
     &		+0.576d0*dsind( U )

	endif

	rLng = DLO+(SLO+DL)/3600.d0
	rLng = mod( mod(rLng,360.d0)+360.d0 , 360.d0 )

	rDis = 10**rDis

	return
	end