       program zody
c        subroutine zody(iracbd, month, iday, iyear, ahour,
c     $     aminute, second, alam, ra, dec, ainu)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C   Calling Program for Ned Wright's Zodiacal Light Routine from
C    Gorjian, Wright, & Chary 2000, ApJ, 536, 550 and
C    Wright 1998, ApJ, 496, 1.  
C   Original program was modifed to account for 
C          1)SIRTF centered coordinates
C          2)IRAC band wavelengths
C          3)IRAC band filters
C   Written by Mark Wolfire
C
C   Input Parameters:
C   ================
C     iracbd  =   IRAC band number (integer) 
C                 0 : program returns intensity at input wavelength 
C                 1 : 3.6 um band
C                 2 : 4.5 um band
C                 3 : 5.6 um band
C                 4 : 8.0 um band
C     month   =  month of observation (integer)
C     iday    =  day of observation  (integer)
C     iyear   =  year of observation (integer)
C     ahour   =  UT hour of observation (real)
C     aminute =  UT minute of observation (real)
C     second  =  UT second of observation (real)
C     ra      =  RA (2000) of observation  (real*4) (degrees)
C     dec     =  Dec (2000) of observation (real*4) (degrees)
C     alam    =  wavelength in microns of observation for iracbd = 0
C                alam is ignored if iracbd equal to 1,2,3,4
C
C
C   Output:
C   ======
C     ainu    = intensity of zodiacal light in MJy sr^-1
C               For iracbd = 0, output is monochromatic intensity
C                 at wavelength alam with no band averaging
C               For iracbd = 1-4, output is monochromatic intensity
C                 at IRAC band center assuming I_nu x nu = constant across
C                 band
C
C   Called Subroutines:
C   ==================
C   sunpos(juliand,ras,decs,elong,oblt) Returns Solar elongation, RA, DEC
C                                         from Julian day number 
C   xyz(juliandr,xs,ys,zs)              Returns Geocentric coordinates of
C                                         Sun from reduced Julian day number
C   euler(ras,decs,gallos,gallas,1)     Returns coordinate transformations
C   interp(JULIANI,YI,NDAY,3,juliand,YSun) Interpolates SIRTF Geocentric
C                                            coordinates for date of obs
C   norm(e)                               Normalizes vectors
C   IPD                                   Ned Wright's code for Zodiacal
C                                            Light calculation
C   iraccol(scattered,thermal,dirbewvl,aracint,iband) Returns band averaged
C                                                  Zodiacal light intensity
C
C   Notes:
C   =====
C   Communication between Ned Wright and Mark Wolfire on the expected 
C   accuracy of the model:
C  
C  I have been comparing the results of your code with those
C  of the Kelsall et code (results kindly supplied by the COBE/DIRBE 
C  folks at Goddard). At a wavelength of 3.5 microns I get a
C  difference of between 5% and approximately 8%. This is for
C  a solar elongation of 120 deg to 80 deg and Galactic latitude = 0,
C  and longitudes of approx 0 to -60 deg. Is this difference entirely
C  reasonable due to your more stingent constraints on the 
C  COBE 12 um emission?
C   Mark
C
C The difference of 5% to 8% is reasonable given the residual shown in the
C 25 micron map with the Kelsall etal model.  So I expect that the
C difference is dominated by the different levels of assumptions about the
C CIRB that went into the models.

C Switching to the SIRTF orbit and SIRTF filter bands will definitely
C increase the uncertainty because these require interpolation for the
C passbands and extrapolation for the orbit change.  I'd be pleased if the
C final zodi count rates observed with SIRTF agree to better than 30% with
C the model predictions.  Once SIRTF is calibrated on orbit, the model will
C probably give reasonable ratios (~10%) for the zodi from place to place.
C Ned Wright
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        parameter (nday=261)
	real r(3)/0.,1.,0./
	real*4 e(3)
	real*8 T68 /1020211225.52D0/
        real*8 juliani(nday),xi(nday),yi(nday),zi(nday),ri(nday)
        real*8 juliand,juliandr,jultime
        real*8 xsun,ysun,zsun
	real ipd
	real denscale,tstar,starlum
	common /hooks/ denscale,tstar,starlum
	dtr = 45./ATAN(1.)
C-----------------------------------------------------------------------------
C   Open file for SIRTF orbit data -> see bottom of data file for 
C     instructions on generating file. 
C               This file is for Launch Date 01.09.03
C-----------------------------------------------------------------------------
        open (unit=10,status='old',file='sirtfld082603.dat')
c        open (unit=10,status='old',file='sirtfld071502.dat')
        do jday = 1,nday
           read(10,*)juliani(jday)
           read(10,*)xi(jday),yi(jday),zi(jday)
           read(10,*)dum,ri(jday),dum
        enddo
        close (unit=10)
C-----------------------------------------------------------------------------
C The following set of input parameters should yield I_nu = 6.053 MJy/sr
          iracbd = 4
          month = 8
          iday = 22
          iyear = 2006
          ahour = 0.0
          aminute = 0.0
          second = 0.0
          alam = 8.0
          ra = 0.0
          dec = 0.0
C-----------------------------------------------------------------------------
C Calculate Julian day and Time from formula in Explanatory Supplement
C   Astronomical Almanac
C-----------------------------------------------------------------------------
          juliand = float((1461*(iyear + 4800 +(month - 14)/12))/4 +
     $        (367*((month-2)-12*((month-14)/12)))/12 -
     $        (3*((iyear+4900 +(month-14)/12)/100))/4+iday-32075)
          jultime = (ahour/24.0d0 - 0.5d0) + (aminute/1440.0d0) + 
     $               (second/86400.0d0)
          juliand = juliand + jultime
          juliandr = juliand - 2400000.0d0
C-----------------------------------------------------------------------------
C Sun and Earth positions 
C-----------------------------------------------------------------------------
          call sunpos(juliand,ras,decs,elong,oblt)
          call xyz(juliandr,xs,ys,zs)
          Rearth = sqrt(xs**2 + ys**2 + zs**2)
          call euler(ras,decs,gallos,gallas,1)
          call euler(ras,decs,ecliptlon,ecliptla,3)
          call euler(ecliptlon,ecliptla,gallos,gallas,5)
          R(1) = cos((elong-180.0)/dtr)*Rearth
          R(2) = sin((elong-180.0)/dtr)*Rearth
          R(3) = 0.0
c          write(6,*)'SUN at l = ',gallos,' b = ',gallas
c          write(6,*)'SUN at RA = ',ras,' DEC = ',decs
c          write(6,*)'SUN at lambda = ',ecliptlon,' beta= ',
c     $        ecliptla
c          write(6,*)'SUN at l = ',gallos,' b = ',gallas
c          write(6,*)'SUN elong = ',elong
c          write(6,*)'Earth Vectors X,Y,Z = ',R(1),R(2),R(3)
C------------------------------------------------------------------------------
C  End of Sun and Earth positions
C------------------------------------------------------------------------------
C  Check if input Julian day is within SIRTF orbit lifetime. If yes
C  use SIRTF centered coordinates. If no use Earth centered coordinates 
C-----------------------------------------------------------------------------
          if(juliand.ge.juliani(1).and.juliand.le.juliani(nday))then
             call interp(juliani,xi,nday,3,juliand,xsun)
             call interp(juliani,yi,nday,3,juliand,ysun)
             call interp(juliani,zi,nday,3,juliand,zsun)
             R(1) = -xsun
             R(2) = -ysun
             R(3) = -zsun
          else
             write(6,*)'Julian day out of range: zodi subtraction'
             write(6,*)'Input Julian day =',juliand
             write(6,*)'First Julian day = ',juliani(1) 
             write(6,*)'Last  Julian day = ',juliani(nday) 
             write(6,*)'Using Earth centered coordinates'          
          end if 
c          write(6,*)'Position Vectors X,Y,Z = ',R(1),R(2),R(3)
          T68 = (juliand-2440000.5D0)*60.0*60.0*24.0
C-----------------------------------------------------------------------------
C  Convert input RA and DEC to ecliptic and Galactic coordinates
C----------------------------------------------------------------------------- 
          call euler(ra,dec,ecliplo,eclipla,3)
          call euler(ra,dec,gallo,galla,1)
C-----------------------------------------------------------------------------
C  Calculate line-of-sight unit vector
C----------------------------------------------------------------------------- 
          E(1) = cos(ecliplo/dtr)*cos(eclipla/dtr)
          E(2) = sin(ecliplo/dtr)*cos(eclipla/dtr)
          E(3) = sin(eclipla/dtr)
          call norm(e) 
          ainu = ipd(E,R,T68,alam,iracbd)
          write(6,*)'ainu',gallo,galla,ra,dec,ainu
	return
	end

	REAL FUNCTION IPD(E,R,T68,LAMBDA,IRACBD)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C   Returns the intensity of the zodiacal light emission in MJy/sr based on
C    Gorjian, Wright, & Chary 2000, ApJ, 536, 550 and
C    Wright 1998, ApJ, 496, 1.  
C
C   Code written by Ned Wright with the following modifications by 
C   Mark Wolfire
C   Original code returned intensity in DIRBE bands to compare directly
C   with observations (i.e. "quoted" values were returned). Intensities
C   at other wavelengths were interpolated between DIRBE band centers.
C   The modifications for SIRTF include 1)setting color corrections 
C   equal to zero (see DIRBE_PLANCK and DIRBE_PLANCK_RANGE) 2)Finding
C   IRAC band averages assuming I_nu x nu = constant. Averages are
C   calculated by interpolating intensity as a function of wavelegnth 
C   between DIRBE band centers.
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	IMPLICIT NONE
	REAL E(3)	! input, unit vector of LOS in ecliptic coordinates
	REAL R(3)	! input, position of observer in ecliptic au's
	REAL*8 T68	! input, time of observation in seconds since 5/24/68
	REAL LAMBDA	! input, microns
C
C	common block for fiddling with emissivities
	REAL ESCAT,EMISS
	COMMON /COEFFS/ ESCAT(10),EMISS(10)
C
	REAL DENSCALE,TSTAR,STARLUM
	COMMON /HOOKS/ DENSCALE,TSTAR,STARLUM
	DATA DENSCALE,TSTAR,STARLUM/1.0,5600.,1.0/
	REAL SAVE_ESCAT(10),SAVE_EMISS(10)
C
C	common block for position
	REAL OBSPOS
	COMMON /POS_OBS/ OBSPOS(3)
C
C       external routines
        REAL PLANCK

c mrm make DUMMY REAL*8 to match with subroutine zmodel
        REAL*8 DUMMY

	REAL SCATTERED(3),THERMAL(8),VAL(10)
	REAL*8 DSCATTERED(3),DTHERMAL(8),DVAL(10)
	INTEGER I
        INTEGER iracbd
        REAL ARACINT
	LOGICAL FIRST /.TRUE./
C
	REAL DIRBEWVL(10)/1.25,2.2,3.5,4.9,12.,25.,60.,100.,140.,240./
	REAL VALUES(8)
        REAL POWER
        REAL SCATVAL,THERMVAL
        REAL TL,TH,TC,RL,RH,RC
C
	IF (FIRST) CALL INIT_ZODI
	FIRST = .FALSE.
	DO I=1,3
	  OBSPOS(I) = R(I)
	ENDDO
	CALL ZMODEL(T68,E,1,10,DVAL,0,DUMMY)
	DO I=1,10
	  VAL(I) = DVAL(I)
	  SAVE_EMISS(I) = EMISS(I)
	  EMISS(I) = 0.
	ENDDO
C	WRITE (*,*) 'VAL = ',VAL
	CALL ZMODEL(T68,E,1,3,DSCATTERED,0,DUMMY)
	DO I=1,3
	  SCATTERED(I) = DSCATTERED(I)
	ENDDO
C	WRITE (*,*) 'SCATTERED=',SCATTERED
	DO I=1,10
	  EMISS(I) = SAVE_EMISS(I)
	  SAVE_ESCAT(I) = ESCAT(I)
	  ESCAT(I) = 0.
	ENDDO
	CALL ZMODEL(T68,E,3,8,DTHERMAL,0,DUMMY)
	DO I=1,8
	  THERMAL(I) = DTHERMAL(I)
	ENDDO
C	WRITE (*,*) 'THERMAL=',THERMAL
	DO I=1,10
	  ESCAT(I) = SAVE_ESCAT(I)
	ENDDO
C-----------------------------------------------------------------------------
C For iracbd = 0 interpolate between DIRBE band centers at wavelength
C    alam. Color corrections = 1.
C For iracbd not equal to 0, return band averaged intensity at band
C    center
C----------------------------------------------------------------------------- 
        if(iracbd.eq.0)then
C
C       divide out by Sun temperature Planck function
C
        DO I=1,3
          SCATTERED(I) = SCATTERED(I)/PLANCK(1.E4/DIRBEWVL(I),5600.)
        ENDDO
C       WRITE (*,*) 'SCATTERED=',SCATTERED
C
C       scale to needed wavelength
        I=2
        IF (LAMBDA.LT.DIRBEWVL(2)) I=1
        POWER = LOG(SCATTERED(I)/SCATTERED(I+1))/
     1       LOG(DIRBEWVL(I)/DIRBEWVL(I+1))
C       WRITE (*,*) I,POWER
        SCATVAL = SCATTERED(2)*PLANCK(1.E4/LAMBDA,5600.)*
     1       (LAMBDA/DIRBEWVL(2))**POWER
C
C       flatten the emissivity
C
        DO I=1,8
          THERMAL(I) = THERMAL(I)*(DIRBEWVL(I+2)+38.)**(2./3.)
        ENDDO
C       WRITE (*,*) 'THERMAL=',THERMAL
        DO I=1,6
          IF (LAMBDA.LT.DIRBEWVL(I+3)) GO TO 100
        ENDDO
C
C       find color temperature of thermal
C
100     TL = 150.
C       WRITE (*,*) I
      RL = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TL)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TL)))
        TH = 350.
      RH = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TH)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TH)))
200     TC = TL*EXP(ALOG(TH/TL)*RL/(RL-RH))
      RC = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TC)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TC)))
        IF (ABS(RC).GT.0.0001) THEN
          IF (ABS(RL).GT.ABS(RH)) THEN
            RL = RC
            TL = TC
          ELSE
            RH = RC
            TH = TC
          ENDIF
          GO TO 200
        ENDIF
C       WRITE (*,*) 'TC,POWER=',TC,POWER
        THERMVAL = (THERMAL(I)/(LAMBDA+38.)**(2./3.))
     1       *PLANCK(1.E4/LAMBDA,TC)/PLANCK(1.E4/DIRBEWVL(I+2),TC)
C       WRITE (*,*) 'SCATVAL,THERMVAL=',SCATVAL,THERMVAL
        IPD = SCATVAL+THERMVAL
C-----------------------------------------------------------------------------
C End of iracbd = 0 
C----------------------------------------------------------------------------
        else
           call iraccol(scattered,thermal,dirbewvl,aracint,iracbd)
           IPD = aracint
        end if
	RETURN
	END

C	"offset widened ringed Good" model with T at 1 au fixed at 280 K
C	breaks from "Gaussian" center to "exponential" tail at sin i[break]
C	with leading/trailing asymmetry and Dermott ring
C	version for bands 1-10
C	CALL INIT_ZODI once before using ZMODEL

	SUBROUTINE INIT_ZODI
C
C	put all the non-linear parameter first
C
C	INCLUDE 'ZMODEL.INC'
	INTEGER NPRM, NLIN, NNLR
	PARAMETER (NPRM=37, NNLR=26, NLIN=NPRM-NNLR)
	INTEGER MINBAND, MAXBAND
	PARAMETER (MINBAND=1, MAXBAND=10)
C
	REAL DENSCALE,TSTAR,STARLUM
	COMMON /HOOKS/ DENSCALE,TSTAR,STARLUM
	COMMON /PARAM/ P(NPRM)
	COMMON /COEFFS/ ESCAT(10),EMISS(10)
C	old values from PHYS
C	REAL SCALE(NLIN)/0.6951908,0.8468719,0.6702332,2.2160680,
C	1	1.4092565,1.2412344,1.2946154,0.8784364,0.7627836,
C	2	0.5358517,0.3957015/
C	REAL P0(NNLR)/1.195541,0.461347,3.701675,0.936077,-1.298938,
C	1	0.371012,-0.083638,-0.042131,0.006141,0.744714,
C	2	5.575954,1.187861,0.001750,0.461642,18.795940,
C	3	-0.033676,0.031110,-0.046420,-0.209133,-0.041613,
C	4	-0.328129,0.590286/

C	values set from fit to COLLECT2.SORT 19-Sep-96:
C	REAL SCALE(NLIN)/0.7045757,0.8581619,0.6577679,2.2966883,
C	1	1.4366271,1.2645168,1.3284225,0.9054061,0.7918135,
C	2	0.5533661,0.4099465/
C	REAL P0(NNLR)/1.200306,0.452485,3.687536,0.924126,-1.296816,
C	1	0.370158,-0.080859,-0.039112,0.005075,0.701273,5.554529,
C	2	1.168095,-0.002007,0.446689,23.744211,-0.034262,0.034704,
C	3	-0.048286,-0.207464,-0.050157,-0.314771,0.574199/
C	values from ZMODEL.FOR_FIZZ3R3, P0 updated to noon 19-Dec-97
C	converged to bands 1-8,VSNZ, 4-Feb-98
C	extended to bands 9 & 10 17-Sep-99
	REAL SCALE(NLIN)/0.6647442,0.7762410,0.7499292,1.8246785,
	1	1.3150289,1.1088078,1.2524812,0.8657632,0.7546864
	2	,0.8147057,0.6130462
	3	/
	REAL P0(NNLR)/1.234614,0.424588,3.578471,0.945027,-1.355945,
	1	0.383785,-0.075830,-0.019477,-0.047109,0.609775,5.540470,
	2	1.385149,0.316134,7.885184,-0.022580,0.028943,-0.026192,
	3	-0.197717,-0.029408,-0.261613,-0.286488,-0.263392,
	4	0.676054,0.531817,0.398778,0.572245/
	IF (NPRM.NE.37) STOP 'BAD NPRM'
	IF (NNLR.NE.26) STOP 'BAD NNLR'
C	WRITE (*,*) 'FIZZ3NP'
	DO I=1,NNLR
	  P(I) = P0(I)
	ENDDO
	DO I=NNLR+1,NPRM
	  P(I) = 1.0
	ENDDO
	DO I=1,10
	  EMISS(I) = 1.E-7
	ENDDO
	DO I=1,10
	  ESCAT(I) = 18.E-14*DIRBE_PLANCK(I,TSTAR)*STARLUM*(5600./TSTAR)**4
	ENDDO
	J = 0
	DO I=MINBAND,3
	  J = J+1
	  ESCAT(I) = SCALE(J)*ESCAT(I)
	ENDDO
	IL = MAX(3,MINBAND)
	DO I=IL,MAXBAND
	  J = J+1
	  EMISS(I) = SCALE(J)*EMISS(I)
	ENDDO
	RETURN
	END

	SUBROUTINE ZMODEL(T68,E,LB,NB,ZOD,DRV,DZ)
C	INCLUDE 'ZMODEL.INC'
	INTEGER NPRM, NLIN, NNLR
	PARAMETER (NPRM=37, NNLR=26, NLIN=NPRM-NNLR)
	INTEGER MINBAND, MAXBAND
	PARAMETER (MINBAND=1, MAXBAND=10)
C
	REAL*8 T68	! time in seconds since 5/24/68 = JD2440000.5
	REAL*4 E(3)	! field of view in ecliptic J2000 unit vector
	INTEGER LB	! first DIRBE band, range 1..10
	INTEGER NB	! number of bands to do
	INTEGER DRV	! 0 for no deriv, 1 for linear only, 2 for all
	REAL*8 ZOD(NB),DZ(NB,NPRM)	! output: zodi in MJy/sr and gradient
C
C	common block for position
	REAL OBSPOS
	COMMON /POS_OBS/ OBSPOS(3)
C
	REAL DENSCALE,TSTAR,STARLUM
	COMMON /HOOKS/ DENSCALE,TSTAR,STARLUM
C
	COMMON /PARAM/ P(NPRM)
	COMMON /COEFFS/ ESCAT(10),EMISS(10)
	REAL*8 SCAT(10),THERM(10),SCAT0(10),THERM0(10)
	REAL*4 R(3),B(3),X(3),POLE(3),OFFSET(3)
	REAL*4 POLD, DELTAP, WATE
C	definitions for longitude of mean Sun
	REAL*8 T0
	PARAMETER (T0 = 682041618.52)
	REAL*4 LONGZ,LONGSLOPE,MEANLONG
	PARAMETER (LONGZ = -77.31159, LONGSLOPE = 359.9562)
C	
	REAL TMAX
	PARAMETER (TMAX=1000.)
C
	INTEGER ISW, ISWMAX
	REAL*4 PLNCK(10),PLISD(10)
	INTEGER NSEG
	INTEGER NPRINT /0/
	REAL*4 WSEG(10)
	REAL*4 FTHERM(10),FSCAT(10)
	REAL*4 LAMBAND(10)/1.25,2.2,3.5,4.9,12.,25.,60.,100.,140.,240./
	REAL PRING(8)/1.705830,0.196297,0.272231,0.251535,0.116663,
	1	-0.269102,-0.932111,1.916413/
C	parameter definitions
	PARAMETER (ND=1)		! density exponent
	PARAMETER (NT=2)		! temperature exponent
	PARAMETER (NDP = 3)		! vertical height scale
	PARAMETER (NDX = 4)		! vertical power law exponent
	PARAMETER (NDW = 5)		! ln (SINIBRK)
	PARAMETER (IX1 = 6, IY1 = 7)	! 10*(zx,zy)
	PARAMETER (IOX1 = 8, IOY1 = 9)	! 10*offset vector
	PARAMETER (IASYM = 10)		! 10*(leading-trailing)
	PARAMETER (IRING = 10)		! Dermott's ring
	PARAMETER (IBAND = 1)		! start of band parameters
	PARAMETER (IPH1 = 20, IPH2 = 23)! phase function coefficients, 2 per band
	PARAMETER (IALB = 26)		! ring relative albedo
	PARAMETER (IT1 = 11)		! temperature for bands at r=1
	PARAMETER (IAMP = 12)		! ring amplitude scaling
C	PARAMETER (IOC = 13)		! parameter for outer cutoff
	PARAMETER (IDMP = 13)		! the "dimple"
	PARAMETER (IDV = 14)		! vertical scale for Dermott
	PARAMETER (NSPH = 15, NSINI2 = 16, 
	1	NG4 = 17, NG16 = 18, NG36 = 19)
	REAL PH1(3),PH2(3)
	EQUIVALENCE (P(IPH1),PH1),(P(IPH2),PH2)
	SAVE
C	get heliocentric position of observer
	DO I=1,3
	  R(I) = OBSPOS(I)
	ENDDO
C	longitude of mean Sun in radians
	MEANLONG = (LONGSLOPE*(T68-T0)/31557600.D0+LONGZ)/57.2957795
	CML = -COS(MEANLONG)
	SML = -SIN(MEANLONG)
c        write(6,*)'mean',meanlong*180.0/3.1415926
c        pause
C	get impact parameter and vector to closest approach of LOS to Sun
	XY = DOT(E,R)		! comp of R parallel to E
	DO I=1,3
	  B(I) = R(I)-XY*E(I)	! comp of R perpendicular to E
	ENDDO
	BB = SQRT(DOT(B,B))
	T1 = XY/BB		! cot
	S = SQRT(1.+T1**2)	! cos = 1/csc
	W1 = 1./(S*(T1+S))
C
C	INTEGRATION VARIABLE IS w = 1 - sin(theta)
C		where X = B + b tan(theta) Y
C		Range of w is from W1 to 0.
C		radius = infinity at w=0
C
	NW=8
	DW=W1/NW
C	offset between emissivity parameter and band number
	KBOFF = 2-MINBAND
	IF (MINBAND.GE.4) KBOFF = KBOFF-1
C	zero the derivatives
	IF (DRV.GT.0) THEN
	  DO IB=1,NB
	    DO JP=1,NPRM
	      DZ(IB,JP) = 0.
	    ENDDO
	  ENDDO
	ENDIF
C	do NNLR+1 evaluations to get value and gradients if requested
	JPMAX = 0
	IF (DRV.GE.2) JPMAX = NNLR
	DO JP=0,JPMAX
	  WATE = 1.0
	  ISWMAX = 1
	  IF (JP.GT.0) THEN
	    DELTAP = 0.015625
	    POLD = P(JP)
	    P(JP) = P(JP)+DELTAP
	    WATE = 0.5/DELTAP
	    ISWMAX = 2
	  ENDIF
	  DO IB=1,NB
	    SCAT(IB) = 0.
	    THERM(IB) = 0.
	  ENDDO
C	computations that don't depend on position
	  POLE(1) = 0.1*P(IX1)
	  POLE(2) = 0.1*P(IY1)
	  POLE(3) = 1.
	  CALL NORM(POLE)
	  OFFSET(1) = 0.1*P(IOX1)
	  OFFSET(2) = 0.1*P(IOY1)
	  OFFSET(3) = 0.
	  SINIBRK = EXP(P(NDW))
	  TEMP1 = EXP(P(IT1))
	  DO ISW=1,ISWMAX
C
C	find the discontinuities in the density model
C	NONE for this model
C
	    NSEG = 1
	    WSEG(1) = 0.
	    WSEG(2) = W1
	    DO ISEG=1,NSEG
	    DW=(WSEG(ISEG+1)-WSEG(ISEG))/NW
	      DO IW=1,NW
	        W = WSEG(ISEG)+DW*(IW-0.5)
	        C=SQRT(W*(2.-W))		!cos(theta)
	        S=1.-W
	        T=BB*S/C
	        DS=BB*DW/C**3
C	FIND THE VECTOR POSITION
	        DO  I=1,3
	          X(I)=B(I)+T*E(I)
	        ENDDO
	        RR=SQRT(DOT(X,X))
C	get position in x,y frame rotating with Earth
	        XD = X(1)*CML+X(2)*SML
	        YD = X(2)*CML-X(1)*SML
C	Dermott's asymmetric cloud
	        DL = ABS(ATAN2(YD,XD)+0.25)
	        IF (DL.GT.0.75) THEN
	          ASYMM = 0.
	        ELSE IF (DL.GT.0.375) THEN
	          ASYMM = (COS(8.37758*DL)-1.)/2 !changed from 8.7758 12-Sep-96
	        ELSE
	          ASYMM = COS(8.37758*DL)	! pi/0.375
	        ENDIF
C	Dermott's ring, with dimple at trailing cloud
	        DERMOTT = -56.5*(SQRT(XD*XD+YD*YD)-1.133
	1		+0.113*P(IDMP)*EXP(-4.*DL**2))**2
	        IF (DERMOTT.LT.-70.) THEN
	          DERMOTT = 0.
	        ELSE
	          DERMOTT = EXP(MIN(10.,DERMOTT))
	        ENDIF
	        SINI = X(3)/RR		! added 14-Sep-96
	        DERMOTT = DERMOTT*EXP(-P(IDV)*SINI**2)
	        ASYMM = DERMOTT*ASYMM
C	PHASE is +1 when looking straight away from Sun, -1 when looking
C	at the Sun, so PHASE is -cosine(scattering angle)
	        PHASE = DOT(E,X)/RR
C	offset: this approximates the semi-major axis.  In the direction of
C	the perihelion, a is greater than R
		RC = RR + DOT(OFFSET,X)
C	tilt
	        SINI = ABS(DOT(X,POLE))/RR
C	note that with a true Good radius, this is really tan|i|
	        IF (SINI.LT.SINIBRK) THEN
	          ZR = (0.5*SINI**2)/SINIBRK
	        ELSE
	          ZR = SINI-0.5*SINIBRK
	        ENDIF
C	temperature depends on RR
	        TEMP = 280.0*EXP(-P(NT)*ALOG(RR/SQRT(STARLUM)))
		ZFACTOR = EXP(-P(NDP)*ZR**P(NDX))
	1		+ P(NSPH) + P(NSINI2)*SINI**2
	2		+ P(NG4)*4.*SINI**2*EXP(-4.*SINI**2)
	3		+ P(NG16)*16.*SINI**2*EXP(-16.*SINI**2)
	4		+ P(NG36)*36.*SINI**2*EXP(-36.*SINI**2)
	        DENSE = (RR/RC)*ZFACTOR*EXP(-P(ND)*ALOG(RC))
	1		* (1.+0.1*(P(IASYM)*ASYMM+P(IRING)*DERMOTT))
C	2		* EXP(MIN(10.,-P(IOC)*RC**2))
		IF (TEMP.GT.TMAX) DENSE=0.
		DENSE = DENSE*DENSCALE
	        SOLAR = DENSE*DS/RR**2
	        THERMAL = DS*DENSE
	        CALL DIRBE_PLANCK_RANGE(3,8,TEMP,PLNCK)
C	        TISD = TEMP	!	400.0/RR**0.4
C	        CALL DIRBE_PLANCK_RANGE(3,8,TISD,PLISD)
C	        DISD = 0.05*P(ISD)
	        DO IB=1,NB
	          JB = LB+IB-1
	          IF (JB.LE.3 .AND. JB.GE.MINBAND) THEN
	             KB = JB-MINBAND+1
	             SCAT(IB) = 
	1		SCAT(IB)+P(KB+NNLR)*ESCAT(JB)*SOLAR
	2		*EXP(PH1(JB)*PHASE+PH2(JB)*PHASE**2)
	          ENDIF
	          IF (JB.GE.3 .AND. JB.LE.MAXBAND) THEN
	            KB = JB+KBOFF
	            THERM(IB)=THERM(IB)+P(KB+NNLR)*EMISS(JB)*
	1		(THERMAL*PLNCK(JB-2))
C	1		(THERMAL*PLNCK(JB-2) +
C	2		DISD*DS*PLISD(JB-2))
	          ENDIF
	        ENDDO
	      ENDDO	! IW
	    ENDDO		! ISEG
	    CALL RING_DING(R,E,PRING(IBAND),PH1,PH2,
	1	TEMP1,FSCAT,FTHERM)
	    DO IB=1,NB
	      JB = LB+IB-1
	      IF (JB.LE.3 .AND. JB.GE.MINBAND) THEN
	        KB = JB-MINBAND+1
                SCAT(IB) = SCAT(IB) + 
	1	ESCAT(JB)*P(KB+NNLR)*P(IALB)*FSCAT(JB)
	      ENDIF
	      IF(JB.GE.3 .AND. JB.LE.MAXBAND) THEN
	        KB = JB+KBOFF
	        THERM(IB) = THERM(IB) + 
	1	EMISS(JB)*P(KB+NNLR)*P(IAMP)*FTHERM(JB)
	      ENDIF
	    ENDDO
	    IF (JP.EQ.0) THEN
	      DO IB=1,NB
	        ZOD(IB) = SCAT(IB)+THERM(IB)
	        SCAT0(IB) = SCAT(IB)
	        THERM0(IB) = THERM(IB)
	      ENDDO
	    ELSE
	      DO IB=1,NB
	        IF (ISW.EQ.1) THEN
	          DZ(IB,JP) = WATE*(SCAT(IB)+THERM(IB))
	        ELSE
	          DZ(IB,JP) = DZ(IB,JP)-WATE*(SCAT(IB)+THERM(IB))
	        ENDIF
	      ENDDO
	      P(JP) = POLD-2.*DELTAP
	    ENDIF
	  ENDDO	! ISW=1...
	IF (JP.GE.1) P(JP) = POLD
	ENDDO	! JP=0,...
C	compute the derivatives wrt the linear parameters
	IF (DRV.GE.1) THEN
	  DO IB=1,NB
	    JB = LB+IB-1
	    IF (JB.LE.3 .AND. JB.GE.MINBAND) THEN
	      KB = JB-MINBAND+1
	      DZ(IB,JB+NNLR) = SCAT0(IB)/P(KB+NNLR)
	    ENDIF
	    IF (JB.GE.3 .AND. JB.LE.MAXBAND) THEN
	      KB = JB+KBOFF
	      DZ(IB,KB+NNLR) = THERM0(IB)/P(KB+NNLR)
	    ENDIF
	  ENDDO
	ENDIF
	RETURN
	END

	SUBROUTINE RING_DING(RIN,EIN,P,PH1,PH2,TEMP1,FSCAT,FTHERM)
C
C	computes the flux in DIRBE bands 1-10 from the zodiacal bands.
C	to be used internally by ZMODEL, so the EMISS and ESCAT scalings
C	are applied in ZMODEL
C
C	RIN is the position of the observer in Sun centered ecliptic
C	coordinates in astronomical units
C	EIN is the unit vector line of sight in ecliptic coordinates
C
C	P is a vector of 8 parameters specific to the rings:
C	P(1) = 10*sin|beta_max| for ring 1, P(2) = 10*optical depth of ring
C	1 material extrapolated to beta=0, r=1
C	P(3-4) are the same for ring 2
C	P(5-6) are 10*Nx/Nz and 10*Ny/Nz for the unit vector N that gives the
C	pole of the ring symmetry plane
C	P(7-8) are 10*Ex and 10*Ey, where (Ex,Ey,Ez=0) is the vector
C	pointing toward the longitude of the perihelion of the forced
C	eccentricity, with magnitude = to the forced eccentricity
C
C	PH1 and PH2 are the phase function parameters for scattering
C	TEMP1 is the temperature at 1 au, assumed to scale like 1/sqrt(r)
C
C	POLE is the pole around which the orbits precess
C	OFFSET is a vector toward the perihelion of the forced
C	eccentricity with a magnitude equal to the foreced
C	eccentricity.
C
C
	REAL DENSCALE,TSTAR,STARLUM
	COMMON /HOOKS/ DENSCALE,TSTAR,STARLUM
C
	REAL*4 R(3),E(3),RIN(3),EIN(3),P(8)
	INTEGER LB,NB
	REAL*4 FSCAT(10),FTHERM(10),PLNCK(10)
	REAL*4 POLE(3),OFFSET(3)
	REAL*4 RMO(3),B(3),X(3)
	REAL*4 STAB(10)
	REAL*4 OMEGA,INC,DTR
	PARAMETER (DTR = 57.29577951)
	REAL PH1(3),PH2(3)
	INTEGER NW
	PARAMETER (NW = 8)
	REAL*4 RRTAB(NW),RCTAB(NW),SINAITAB(NW)
	REAL*4 AB(2)/		0.03,		0.03/	! amplitudes
	REAL*4 SINI(2)/		0.0247811,	0.176741/ ! sin(inc)
	REAL*4 RMAX(2,2)/	3.14,		3.02, 	! outer radii
	1			0.1,		0.1/	! inner radii
	DATA IFENT /0/
	SECI = SQRT(1.+0.01*(P(5)**2+P(6)**2))
	POLE(1) = 0.1*P(5)/SECI
	POLE(2) = 0.1*P(6)/SECI
	POLE(3) = 1./SECI
	OFFSET(1) = 0.1*P(7)
	OFFSET(2) = 0.1*P(8)
	OFFSET(3) = 0.0
	SINI(1) = 0.1*P(1)
	AB(1) = 0.1*P(2)
	SINI(2) = 0.1*P(3)
	AB(2) = 0.1*P(4)
	LB = 1
	NB = 10
	DO IB=1,NB
	  FSCAT(IB)=0.
	  FTHERM(IB)=0.
	ENDDO
	DO I=1,3
	  E(I) = EIN(I)
	  R(I) = RIN(I)
	ENDDO
	CALL NORM(E)
C	get impact parameter and vector to closest approach of LOS to Sun
	XY = DOT(E,R)		! comp of R parallel to E
	DO I=1,3
	  B(I) = R(I)-XY*E(I)	! comp of R perpendicular to E
	ENDDO
	BB2 = DOT(B,B)
	BB = SQRT(BB2)
	T1 = XY/BB		! cot
C	R = B + BB*T1*E
	S = SQRT(1.+T1**2)	! cos = 1/csc
	W1 = 1./(S*(T1+S))
	DO I=1,3
	  RMO(I) = B(I)-OFFSET(I)
	ENDDO
	BP = DOT(B,POLE)
	BO = DOT(B,OFFSET)
	EO = DOT(E,OFFSET)
	EP = DOT(E,POLE)
C	Equation for intersection of LOS and cone
C	(B+E*s)*Pole = z = +/-|sin i|*|B+E*s|
C	[(E*Pole)s + B*Pole]^2 = (sin i)^2*|B+E*s|^2
C	[(E*Pole)^2 - (sin i)^2]s^2 
C	+2*[(E*Pole)*{B*Pole} - (sin i)^2*{E*B}]s
C	+[{B*Pole}^2 - (sin i)^2*B^2] = 0
	DO IR=1,2
	  NSEG = 0
	  A = EP**2-SINI(IR)**2
	  BQ = 2.0*EP*BP
	  C = BP**2-BB2*SINI(IR)**2
	  DISC = BQ*BQ-4.*A*C
	  IF (DISC.GT.0.) THEN
	    DISC = SQRT(DISC)
	    S = (-BQ+DISC)/(2.*A)
	    IF (S.GT.XY) THEN
	      NSEG = NSEG+1
	      STAB(NSEG) = S
	    ENDIF
	    S = (-BQ-DISC)/(2.*A)
	    IF (S.GT.XY) THEN
	      NSEG = NSEG+1
	      STAB(NSEG) = S
	    ENDIF
	  ELSE IF (DISC.EQ.0.) THEN
	    S = -BQ/(2.*A)
	    IF (S.GT.XY) THEN
	      NSEG = NSEG+1
	      STAB(NSEG) = S
	    ENDIF
	  ENDIF
C	Equation for intersection of LOS and surface of constant semi-major
C	axis: a = |R| + Offset*R = RMAX
C	(RMAX - Offset*R)^2 = |R|^2
C	RMAX^2 - 2*RMAX*(Offset*B + s*(Offset*E)) 
C	+ (Offset*B)^2 + 2*s*(Offset*B)*(Offset*E) + s^2* (Offset*E)^2
C	= (B+E*s)^2 = B*B + s^2
C
C	(1-(Offset*E)^2)*s^2 +
C	(2*RMAX*(Offset*E)-2*(Offset*B)*(Offset*E))*s+
C	B*B - (RMAX-(Offset*B))^2 = 0
C
	  DO INOUT=1,2
	    A = 1.-EO**2
	    BQ = 2.*(RMAX(IR,INOUT)-BO)*EO
	    C = BB2-(RMAX(IR,INOUT)-BO)**2
	    DISC = BQ*BQ-4.*A*C
	    IF (DISC.GT.0.) THEN
	      DISC = SQRT(DISC)
	      S = (-BQ+DISC)/(2.*A)
	      IF (S.GT.XY) THEN
	        NSEG = NSEG+1
	        STAB(NSEG) = S
	      ENDIF
	      S = (-BQ-DISC)/(2.*A)
	      IF (S.GT.XY) THEN
	        NSEG = NSEG+1
	        STAB(NSEG) = S
	      ENDIF
	    ELSE IF (DISC.EQ.0.) THEN
	      S = -BQ/(2.*A)
	      IF (S.GT.XY) THEN
	        NSEG = NSEG+1
	        STAB(NSEG) = S
	      ENDIF
	    ENDIF
	  ENDDO	! INOUT=1,2
C	sort the S table into ascending order, note that NSEG <=6
	  DO IS=1,(NSEG-1)
	    DO JS=(IS+1),NSEG
	      IF (STAB(IS).GT.STAB(JS)) THEN
	        W = STAB(IS)
	        STAB(IS) = STAB(JS)
	        STAB(JS) = W
	      ENDIF
	    ENDDO
	  ENDDO
C	do the integrals
C
C	INTEGRATION VARIABLE IS w = 1 - sin(theta)
C		where X = B + bb tan(theta) E
C		Range of w is from W1 to 0.
C		radius = infinity at w=0
C	note that STAB has values of s = b tan(theta) where discontinuities
C	occur.  w = 1 - tan(theta)/sec(theta) = 1-(s/b)/sqrt(1+(s/b)^2)
C
	  STAB(NSEG+2) = 0.
	  DO JS=1,NSEG
	    IS = NSEG+1-JS
	    STAB(IS+1) = 1.-(STAB(IS)/BB)/SQRT(1.+(STAB(IS)/BB)**2)
	  ENDDO
	  NSEG = NSEG+1
	  STAB(1) = W1
	  DO IS=1,NSEG
	    DW=(STAB(IS)-STAB(IS+1))/NW
	    IF (DW.LE.0.) GO TO 90
	    ICNT = 0
	    DO IW=1,NW
	      W = STAB(IS)-DW*(IW-0.5)
	      C = W*(2.-W)
	      IF (C.LE.0.0) THEN
C	        WRITE (*,*) 'W =',W,IW,IS
C	        WRITE (*,FMT='(8F10.6)') (STAB(III+1),III=0,NSEG)
	        GO TO 80
	      ENDIF
	      C=SQRT(C)			! cos(theta)
	      S=1.-W
	      T=BB*S/C
	      DS=BB*DW/C**3
C	FIND THE VECTOR POSITION
	      DO I=1,3
	        X(I)=B(I)+T*E(I)
	      ENDDO
	      RR=SQRT(DOT(X,X))
C	PHASE is +1 when looking straight away from Sun, -1 when looking
C	at the Sun, so PHASE is -cosine(scattering angle)
	      PHASE = DOT(E,X)/RR
C	temperature depends on RR
	      TEMP = TEMP1/SQRT(RR/SQRT(STARLUM))
C	use offset, tilted cloud for density
	      Z = ABS(DOT(X,POLE))
	      SINAI = Z/RR
C	approximate semi-major axis is RC
	      RC = RR + DOT(OFFSET,X)
	      RRTAB(IW) = RR
	      RCTAB(IW) = RC
	      SINAITAB(IW) = SINAI
	      IF (RC.GT.RMAX(IR,1) .OR. RC.LT.RMAX(IR,2)
	1	.OR. SINAI.GT.SINI(IR)) THEN
	        GO TO 90	! bail out to save time
	      ELSE
C	density factor (r/a)/a^alpha = RR/RC**2
	      DENSE = AB(IR)*COSH(1.719266*SINAI/SINI(IR))*(RR/RC**2)
	      DENSE = DENSCALE*DENSE
	      THERMAL = DS*DENSE
	      CALL DIRBE_PLANCK_RANGE(LB,NB,TEMP,PLNCK)
	      SOLAR = DENSE*DS/RR**2
	      DO IB=1,NB
	        JB = LB+IB-1
	        FSCAT(IB) = FSCAT(IB)+SOLAR
	1	* EXP(PH1(MIN(3,IB))*PHASE+PH2(MIN(3,IB))*PHASE**2)
	        FTHERM(IB) = FTHERM(IB)+THERMAL*PLNCK(IB)
	      ENDDO
	      ENDIF
80	      CONTINUE
	    ENDDO	! IW=1,NW
90	    CONTINUE	! bailout for zero density segments
	  ENDDO		! IS=1,NSEG
	ENDDO	! IR=1,2
	RETURN
	END

	REAL*4 FUNCTION DIRBE_PLANCK(BAND,T)
	REAL DIRBEWVL(10)/1.25,2.2,3.5,4.9,12.,25.,60.,100.,140.,240./

C
C	computes response of DIRBE band "BAND" to a blackbody of temperature
C	T in Kelvin
C	answer is in MJy/sr
C
	INTEGER*4 BAND
	REAL*4 T
	REAL*4 VALUE
C--------------------------------------------------------------------------
C DIRBE_PLANK_RANGE is for K(T)xB(T) for the DIRBE BANDs. MGW.
C   Calling parameters are (1st DIRBE BAND, Number of Bands to Calculate,
C   Temperature, Values)  
C-------------------------------------------------------------------------- 
c	CALL DIRBE_PLANCK_RANGE(BAND,1,T,VALUE)
c	DIRBE_PLANCK = VALUE
C--------------------------------------------------------------------------
C The Following is a substitute routine to find I_nu (MJy/sr)  MGW
C   This needs to be replaced by an IRAC integrated intensity ?
C   Nope - it cannot be an IRAC integrated intensity because parameter
C   values are known only at DIRBE bands.
C-------------------------------------------------------------------------
        DIRBE_PLANCK = PLANCKMJ(1.E4/DIRBEWVL(BAND),T)
	RETURN
	END

	SUBROUTINE DIRBE_PLANCK_RANGE(BAND,NB,T,VALUES)
C
C	computes response of DIRBE band "BAND" to a blackbody of temperature
C	T in Kelvin
C	answer is in MJy/sr
C
	REAL DIRBEWVL(10)/1.25,2.2,3.5,4.9,12.,25.,60.,100.,140.,240./
	INTEGER*4 BAND,NB
	REAL*4 T
	REAL*4 VALUES(NB)
	PARAMETER (T0 = 2.726, NT = 921, NBN = 10)
	REAL*4 RESPONSE(NBN,NT)
	DATA IFENT /0/
	SAVE
C-----------------------------------------------------------------------------
C  DIRBE response [K(T)xB(T)] is returned only if "quoted" DIRBE values 
C  are wanted. We return "true" B(T) and interpolate between DIRBE
C  wavebands to find band average IRAC intensities
C  The following silly statement is included to skip existing code.
C----------------------------------------------------------------------------  
        if(4.gt.2)goto 100
	IF (IFENT.EQ.0) THEN
	  IFENT = 27
	  OPEN (UNIT=27,
	1	FILE='DIRBE_PLANCK_TABLE.DAT',
	2	STATUS='OLD')
	  DO I=1,NT
	    READ (27,FMT='(F11.2,10F10.6)') TT,(RESPONSE(J,I),J=1,NBN)
	    DO J=1,NBN
	      IF (RESPONSE(J,I).LT.-75.) THEN
	        RESPONSE(J,I) = 0.
	      ELSE
	        RESPONSE(J,I) = EXP(RESPONSE(J,I))
	      ENDIF
	    ENDDO
	  ENDDO
	  CLOSE (UNIT=27)
	ENDIF
	IF (BAND.LT.1 .OR. BAND+NB-1.GT.NBN) THEN
	  WRITE (*,*) 'BAD BAND IN DIRBE_PLANCK_RANGE'
	  BAND = SQRT(25.-(BAND-5.5)**2)
	ENDIF
	IF (T.LT.T0) THEN
	  DIRBE_PLANCK = 0.
	  RETURN
	ENDIF
C	note that first element of table is for T0, so when
	F = 230.*ALOG10(T/T0)
C	F = 0+ we should interpolate using rows 0,1,2 and 3
	I = F
	F = F-I
C	coefficients for cubic interpolation
	F0 = -0.16666667*F*(F-1.)*(F-2.)
	F1 = 0.5*(F*F-1.)*(F-2.)
	F2 = -0.5*(F+1.)*F*(F-2.)
	F3 = 0.16666667*F*(F*F-1.)
 100	continue
	DO IB = 1,NB
	  JB = BAND+IB-1
c	  VALUES(IB) = F0*RESPONSE(JB,I) + F1*RESPONSE(JB,I+1)
c	1	+ F2*RESPONSE(JB,I+2) + F3*RESPONSE(JB,I+3)
C---------------------------------------------------------------------------
C Routine to return true I_nu MJy/sr in DIRBE BANDS instead of
C      K(T) x B(T)
C --------------------------------------------------------------------------
          VALUES(IB) = PLANCKMJ(1.E4/DIRBEWVL(JB),T)
	ENDDO
	RETURN
	END

	FUNCTION DOT(A,B)
	DIMENSION A(3),B(3)
	DOT=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
	RETURN
	END

	SUBROUTINE NORM(A)
C
C	NORMALIZES A THREE VECTOR
C
	DIMENSION A(3)
	AM=MAX(ABS(A(1)),ABS(A(2)),ABS(A(3)))
	IF(AM.GT.1.E10) THEN
	  DO I=1,3
	    A(I)=A(I)/AM
	  ENDDO
	ENDIF
	AA=SQRT(A(1)**2+A(2)**2+A(3)**2)
	DO I=1,3
	  A(I)=A(I)/AA
	ENDDO
	RETURN
	END

	SUBROUTINE CROSS(A,B,C)
	DIMENSION A(3),B(3),C(3)
C
C	CROSS PRODUCT C=A X B
C
	C(3)=A(1)*B(2)-A(2)*B(1)
	C(1)=A(2)*B(3)-A(3)*B(2)
	C(2)=A(3)*B(1)-A(1)*B(3)
	RETURN
	END

	SUBROUTINE ETOC(E,C)
	REAL*4 E(3),C(3),TEMP
C
C	E	INPUT	vector in ecliptic coordinates (J2000.0)
C	C	OUTPUT	vector in celestial coordinates (J2000.0)
C
	PARAMETER (COBL=.917482058,SOBL=.39777716)
	C(1)=E(1)
	TEMP=COBL*E(3)+SOBL*E(2)
	C(2)=COBL*E(2)-SOBL*E(3)
	C(3)=TEMP
	RETURN
	END

	SUBROUTINE MAV(A,VSC,VCEL)
C
C	MULTIPLIES MATRIX A TIMES VECTOR VSC GIVING VECTOR VCEL
C
	REAL*4 A(3,3),VSC(3),VCEL(3),SC(3),VC
	DO 5 I=1,3
5	SC(I)=VSC(I)
	DO 15 I=1,3
	VC=0.
	DO 10 J=1,3
10	VC=VC+A(I,J)*SC(J)
15	VCEL(I)=VC
	RETURN
	END

	SUBROUTINE VMAT(A,VCEL,VSC)
C
C	MULTIPLIES TRANSPOSE OF A TIMES VCEL GIVING VSC
C
	REAL*4 A(3,3),VCEL(3),VSC(3),CEL(3),SC
	DO 5 I=1,3
5	CEL(I)=VCEL(I)
	DO 15 I=1,3
	SC=0.
	DO 10 J=1,3
10	SC=SC+A(J,I)*CEL(J)
15	VSC(I)=SC
	RETURN
	END

	FUNCTION PLANCK(ANU,T)
C
C	GIVES ERGS/CM**2/SR/SEC/WAVENUMBER
C	T IN KELVIN
C	ANU IN 1/CM
C
	REAL*4 THCC,HCK
	PARAMETER (THCC=1.1910439E-5,HCK=1.438769)
	IF(T.GT.0.) GO TO 10
5	PLANCK=0.
	RETURN
10	X=HCK*ANU/T
	IF(X.GT.50.) GO TO 5
	IF(X.GT..01) GO TO 15
	D=X*(1.+X*(.5+X/6.))
	GO TO 20
15	D=EXP(X)-1.
20	PLANCK=THCC*ANU**3/D
	RETURN
	END



	FUNCTION PLANCKMJ(ANU,T)
C
C	GIVES MJy/Sr
C	T IN KELVIN
C	ANU IN 1/CM
C
	REAL*4 THCC,HCK,SCALE
	PARAMETER (THCC=1.1910439E-5,HCK=1.438769)
        SCALE = 1.0E-26*1.0e7/1.0e4/ANU*1.0e6*
     1          (2.997925e10*ANU) 
	IF(T.GT.0.) GO TO 10
5	PLANCKMJ=0.
	RETURN
10	X=HCK*ANU/T
	IF(X.GT.50.) GO TO 5
	IF(X.GT..01) GO TO 15
	D=X*(1.+X*(.5+X/6.))
	GO TO 20
15	D=EXP(X)-1.
20	PLANCKMJ=THCC*ANU**3/D/SCALE
	RETURN
	END

      subroutine iraccol(scattered,thermal,dirbewvl,aracint,iband)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C  Subroutine to return IRAC band averaged intensity in MJy/sr
C  written by Mark Wolfire
C
C  Input Parameters:
C  ================
C    scattered = Scattered light component of zodiacal light at first three
C                 DIRBE wavebands (MJy/sr) (1.25, 2.2, 3.5 microns)
C    thermal   = Thermal component of zodiacal light at last eight DIRBE
C                  wavebands (MJy/sr) (3.5,4.9,12,25,60,100,140,240 microns) 
C    dirbewvl  = DIRBE band centers (microns)
C    iband     = IRAC band number
C                1 : 3.6 um band
C                2 : 4.5 um band
C                3 : 5.6 um band
C                4 : 8.0 um band
C
C  Outputs:
C  =======
C    aracint  = intensity at IRAC band center assuming I_nu x nu = constant
C                (MJy/sr) i.e. "quoted" intensities where 
C               I_nu (quoted) = K x I_nu (actual) and K is calculated
C               as in SIRTF Cailibration: Units and Conventions by
C               Bill Reach.  
C               Integration constants = integral [ (nu_0/nu) R/nu dnu] are 
C               precalculated and stored in bandint.
C
C   Notes:
C   =====
C   IRAC band filters R are read in as needed 
C   Assumes filters are electrons/photon
C        irac_tr1.dat = IRAC band 1
C        irac_tr2.dat = IRAC band 2
C        irac_tr3.dat = IRAC band 3
C        irac_tr4.dat = IRAC band 4
C   Intensity is interpolated between DIRBE band centers using algorithm
C     from Ned Wright.
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        dimension scattered(3),thermal(8),dirbewvl(10)
        dimension wave(201),tran(201),aint(201),ainu(201)
        dimension bandint(4)
        character*12 datafile
        character*1 ib
        integer iband
        real lambda,ipd
        integer i 
	REAL POWER
	REAL SCATVAL,THERMVAL
	REAL TL,TH,TC,RL,RH,RC
        
C
C	external routines
	REAL PLANCK

C-----------------------------------------------------------------------------
C   Precalculated integration constants for each IRAC band  (07/18/02)
C   Assumes response function is eletrons/photon and
C   IRAC 1 lambda_0 = 3.5482
C   IRAC 2 lambda_0 = 4.4916
C   IRAC 3 lambda_0 = 5.6603
C   IRAC 4 lambda_0 = 7.8704
C----------------------------------------------------------------------------
C Constants for R
        data bandint/0.126918256,0.143180043,0.064407073,0.139533788/
C Constants for R_nu 
c        data bandint/1.08712E+13,9.53873E12,3.32909E12,5.22890E+12/
        c = 2.997925e14
C----------------------------------------------------------------------------
C   Open IRAC band response data file
C----------------------------------------------------------------------------
        write(ib,'(I1)'),iband
        datafile = 'irac_tr'//ib//'.dat'
        open(unit=20,status='old',file=datafile)
C
C       divide out by Sun temperature Planck function
C
        DO I=1,3
          SCATTERED(I) = SCATTERED(I)/PLANCK(1.E4/DIRBEWVL(I),5600.)
        ENDDO
c       WRITE (6,*) 'SCATTERED=',SCATTERED
C
C       flatten the emissivity
C
        DO I=1,8
          THERMAL(I) = THERMAL(I)*(DIRBEWVL(I+2)+38.)**(2./3.)
        ENDDO
C       WRITE (*,*) 'THERMAL=',THERMAL

C----------------------------------------------------------------------------
C  Start integration loop for each wavelength in filter response file
C  Loop calculates I_nu x R where I_nu is intensity interpolated between
C  DIRBE bands 
C----------------------------------------------------------------------------
        read(20,*)nlines
        do il = 1,nlines
           inu = nlines + 1 -il 
           read(20,*)wave(il),tran(inu)
           ainu(inu) = c/wave(il)
           lambda = wave(il)
C
C       scale to needed wavelength
        I=2
        IF (LAMBDA.LT.DIRBEWVL(2)) I=1
        POWER = LOG(SCATTERED(I)/SCATTERED(I+1))/
     1       LOG(DIRBEWVL(I)/DIRBEWVL(I+1))
c       WRITE (*,*) I,POWER
        SCATVAL = SCATTERED(2)*PLANCK(1.E4/LAMBDA,5600.)*
     1       (LAMBDA/DIRBEWVL(2))**POWER

        DO I=1,6
          IF (LAMBDA.LT.DIRBEWVL(I+3)) GO TO 100
        ENDDO
C
C       find color temperature of thermal
C
100     TL = 150.
C       WRITE (*,*) I
      RL = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TL)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TL)))
      TH = 350.
      RH = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TH)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TH)))
200   TC = TL*EXP(ALOG(TH/TL)*RL/(RL-RH))
      RC = LOG((THERMAL(I)/THERMAL(I+1))/(PLANCK(1.E4/DIRBEWVL(I+2),TC)/
     1       PLANCK(1.E4/DIRBEWVL(I+3),TC)))
       IF (ABS(RC).GT.0.0001) THEN
          IF (ABS(RL).GT.ABS(RH)) THEN
            RL = RC
            TL = TC
          ELSE
            RH = RC
            TH = TC
          ENDIF
          GO TO 200
        ENDIF
C       WRITE (*,*) 'TC,POWER=',TC,POWER
        THERMVAL = (THERMAL(I)/(LAMBDA+38.)**(2./3.))
     1       *PLANCK(1.E4/LAMBDA,TC)/PLANCK(1.E4/DIRBEWVL(I+2),TC)
C       WRITE (*,*) 'SCATVAL,THERMVAL=',SCATVAL,THERMVAL
        IPD = SCATVAL+THERMVAL
C Integral assumes response function is electron/photon
        aint(inu) = IPD*tran(inu)/ainu(inu)
C Integral assumes response function is per Hz
c        aint(inu) = IPD*tran(inu)
        end do
C-----------------------------------------------------------------------------
C  End I_nu x R loop
C-----------------------------------------------------------------------------
        call oneteg(ainu,aint,aracint,nlines)
        aracint = aracint/bandint(iband)
C----------------------------------------------------------------------------
C  aracint = returned intensity
C---------------------------------------------------------------------------- 
        close(20)
        return
        end

	SUBROUTINE ONETEG(X,Y,TEG,NPOINT)
	DIMENSION X(1),Y(1)
	TEG = 0.0
	DO 10 I=2,NPOINT
	TEG = TEG + 0.5*(Y(I)+Y(I-1))*(X(I)-X(I-1))	
10	CONTINUE
	RETURN
	END
      subroutine sunpos( jd, ra, dec, longmed, oblt)
c+
c NAME:
c       SUNPOS
c PURPOSE:
c       To compute the RA and Dec of the Sun at a given date.
c
c CALLING SEQUENCE:
c       SUNPOS, jd, ra, dec, [elong, obliquity, /RADIAN ]
c INPUTS:
c       jd    - The Julian date of the day (and time), scalar or vector
c               usually double precision
c OUTPUTS:
c       ra    - The right ascension of the sun at that date in DEGREES
c               double precision, same number of elements as jd
c       dec   - The declination of the sun at that date in DEGREES
c
c OPTIONAL OUTPUTS:
c       elong - Ecliptic longitude of the sun at that date in DEGREES.
c       obliquity - the obliquity of the ecliptic, in DEGREES
c
c OPTIONAL INPUT KEYWORD:
c       /RADIAN - If this keyword is set and non-zero, then all output variables 
c               are given in Radians rather than Degrees
c
c NOTES:
c       The accuracy in the 20th century  should be within 1"c however this 
c       has not been extensively tested.
c
c       The returned RA and Dec are in the given date's equinox.
c
c       Procedure was extensively revised in May 1996, and the new calling
c       sequence is incompatible with the old one.
c METHOD:
c       Uses a truncated version of Newcomb's Sun.    Adapted from the IDL
c       routine SUN_POS by CD Pike, which was adapted from a FORTRAN routine
c       by B. Emerson (RGO).
c EXAMPLE:
c       (1) Find the apparent RA and Dec of the Sun on May 1, 1982
c       
c       IDL> jdcnv, 1982, 5, 1,0 ,jd      cFind Julian date jd = 2445090.5   
c       IDL> sunpos, jd, ra, dec
c       IDL> print,adstring(ra,dec,2)
c                02 31 32.61  +14 54 34.9
c
c       The Astronomical Almanac gives 02 31 32.58 +14 54 34.9 so the error
c               in SUNPOS for this case is < 0.5".      
c
c       (2) Find the apparent RA and Dec of the Sun for every day in 1997
c
c       IDL> jdcnv, 1997,1,1,0, jd                cJulian date on Jan 1, 1997
c       IDL> sunpos, jd+ dindgen(365), ra, dec    cRA and Dec for each day 
c
c MODIFICATION HISTORY:
c       Written by Michael R. Greason, STX, 28 October 1988.
c       Accept vector arguments, W. Landsman     April,1989
c       Eliminated negative right ascensions.  MRG, Hughes STX, 6 May 1992.
c       Rewritten using the 1993 Almanac.  Keywords added.  MRG, HSTX, 
c               10 February 1994.
c       Major rewrite, improved accuracy, always return values in degrees
c       W. Landsman  May, 1996 
c       Added /RADIAN keyword,    W. Landsman       August, 1997
c       Converted to IDL V5.0   W. Landsman   September 1997
c       Converted to FORTRAN!, by M. Wolfire for use in SIRTF GLIMPSE 
c           pipeline February 2001
c-
c On_error,2
c                       Check arguments.
c if N_params() LT 3 then begin 
c     print, 'Syntax - SUNPOS, jd, ra, dec, [elong, obliquity, /RADIAN] '
c     print, 'Inputs  -  jd (Julian date)'
c     print, 'Outputs - Apparent RA and Dec, longitude, & obliquity'
c     print, 'All angles in DEGREES unless /RADIAN is set'
c     return
c endif
      real*8 jd,t,l,me,mv,mm,mj,marscorr,mooncorr,longterm,jupcorr
      real*4 longmed
      pi = 3.1415927
      dtor = pi/180.0d0       !(degrees to radian, double precision)

c  form time in Julian centuries from 1900.0

      t = (jd - 2415020.0d0)/36525.0d0

c  form sun's mean longitude

      l = (279.696678d0+mod((36000.768925d0*t),360.0d0))*3600.0d0
c  allow for ellipticity of the orbit (equation of centre)
c  using the Earth's mean anomoly ME

      me = 358.475844d0 + mod((35999.049750D0*t),360.0d0)
      ellcor  = (6910.1d0 - 17.2D0*t)*sin(me*dtor) + 
     $           72.3D0*sin(2.0D0*me*dtor)
      l = l + ellcor

c allow for the Venus perturbations using the mean anomaly of Venus MV

      mv = 212.603219d0 + mod((58517.803875d0*t),360.0d0) 
      vencorr = 4.8D0 * cos((299.1017d0 + mv - me)*dtor) + 
     $   5.5D0 * cos((148.3133d0 +  2.0D0 * mv  -  2.0D0 * me )*dtor) +
     $   2.5D0 * cos((315.9433d0 +  2.0D0 * mv  -  3.0D0 * me )*dtor) +
     $   1.6D0 * cos((345.2533d0 +  3.0D0 * mv  -  4.0D0 * me )*dtor) +
     $   1.0D0 * cos((318.15d0   +  3.0D0 * mv  -  5.0D0 * me )*dtor)
      l = l + vencorr

c  Allow for the Mars perturbations using the mean anomaly of Mars MM

      mm = 319.529425d0  +  mod(( 19139.858500d0 * t),360.0d0 )
      marscorr = 2.0d0 * cos((343.8883d0 -  2.0d0 * mm  +  2.0d0* 
     $    me)*dtor)+1.8D0 * cos((200.4017d0 -  2.0d0 * mm  + me) 
     $    * dtor)
      l = l + marscorr

c Allow for the Jupiter perturbations using the mean anomaly of
c Jupiter MJ

      mj = 225.328328d0  +  mod(( 3034.6920239d0 * t),360.0d0 )
      jupcorr = 7.2d0 * cos(( 179.5317d0 - mj + me )*dtor) + 
     $   2.6d0 * cos((263.2167d0  -  MJ ) *dtor) + 
     $   2.7d0 * cos(( 87.1450d0  -  2.0d0 * mj  +  2.0D0 * me ) 
     $   *dtor) + 1.6d0 * cos((109.4933d0  -  2.0d0 * mj  +  me ) 
     $   *dtor)
      l = l + jupcorr
c Allow for the Moons perturbations using the mean elongation of
c the Moon from the Sun D

      d = 350.7376814d0  + mod(( 445267.11422d0 * t),360.0d0 )
      mooncorr  = 6.5d0 * sin(d*dtor)
      l = l + mooncorr
c Allow for long period terms

       longterm  = + 6.4d0 * sin(( 231.19d0  +  20.20d0 * t )*dtor)
       l  =    l + longterm
       l  =  mod(( l + 2592000.0d0),1296000.0d0) 
       longmed = l/3600.0d0
c Allow for Aberration

       l  =  l - 20.5d0

c Allow for Nutation using the longitude of the Moons mean node OMEGA

       omega = 259.183275d0 - mod(( 1934.142008d0 * t ),360.0d0 )
       l  =  l - 17.2d0 * sin(omega*dtor)

c Form the True Obliquity

       oblt  = 23.452294d0 - 0.0130125d0*t + 
     $      (9.2d0*cos(omega*dtor))/3600.0d0

c Form Right Ascension and Declination

       l = l/3600.0d0
       ra  = atan2( sin(l*dtor) * cos(oblt*dtor) , cos(l*dtor) )
c       ra  = atan( sin(l*dtor) * cos(oblt*dtor))
       neg = 0
       if(ra.LT.0.0d0)ra = ra + 2.0d0*pi
c       neg = where(ra LT 0.0d0, Nneg) 
c       if Nneg GT 0 then ra[neg] = ra[neg] + 2.0d*pi

       dec = asin(sin(l*dtor) * sin(oblt*dtor))
 
c if keyword_set(RADIAN) then begin
c        oblt = oblt*dtor 
c        longmed = longmed*dtor
c endif else begin
        ra = ra/dtor
        dec = dec/dtor
c endelse
      return
      end
      SUBROUTINE EULER(AI,BI,AO,BO,ISELECT)
c+
c NAME:
c     EULER
c PURPOSE:
c     Transform between Galactic, celestial, and ecliptic coordinates.
c EXPLANATION:
c     Use the procedure ASTRO to use this routine interactively
c
c CALLING SEQUENCE:
c      EULER, AI, BI, AO, BO, [ SELECT, /FK4 ] 
c
c INPUTS:
c       AI - Input Longitude in DEGREES, scalar or vector.  If only two 
c               parameters are supplied, then  AI and BI will be modified to 
c               contain the output longitude and latitude.
c       BI - Input Latitude in DEGREES
c
c OPTIONAL INPUT:
c       SELECT - Integer (1-6) specifying type of coordinate transformation.  
c
c      SELECT   From          To        |   SELECT      From            To
c       1     RA-Dec (2000)  Galactic   |     4       Ecliptic      RA-Dec    
c       2     Galactic       RA-DEC     |     5       Ecliptic      Galactic  
c       3     RA-Dec         Ecliptic   |     6       Galactic      Ecliptic  
c
c      If omitted, program will prompt for the value of SELECT
c      Celestial coordinates (RA, Dec) should be given in equinox J2000 
c      unless the /FK4 keyword is set.
c OUTPUTS:
c       AO - Output Longitude in DEGREES
c       BO - Output Latitude in DEGREES
c
c INPUT KEYWORD:
c       /FK4 - If this keyword is set and non-zero, then input and output 
c             celestial and ecliptic coordinates should be given in equinox 
c             B1950.
c
c NOTES:
c       EULER was changed in December 1998 to use J2000 coordinates as the 
c       default, ** and may be incompatible with earlier versions***.
c REVISION HISTORY:
c       Written W. Landsman,  February 1987
c       Adapted from Fortran by Daryl Yentis NRL
c       Converted to IDL V5.0   W. Landsman   September 1997
c       Made J2000 the default, added /FK4 keyword  W. Landsman December 1998
c       Converted back to FORTRAN!, by M. Wolfire for use in SIRTF 
c           GLIMPSE pipeline  February 2001
c-
c On_error,2

c npar = N_params()
c if npar LT 4 then begin
c    print,'Syntax - EULER, AI, BI, A0, B0, [ SELECT, /FK4 ]'
c    print,'    AI,BI - Input longitude,latitude in degrees'
c    print,'    AO,BO - Output longitude, latitude in degrees'
c    print,'    SELECT - Scalar (1-6) specifying transformation type'
c    return
c endif
      real*8 psi(6),stheta(6),ctheta(6),phi(6)
c      real*4 ai,ao,bi,bo
      data psi/ 0.57477043300D0, 4.9368292465D0,  
     $       0.00000000000D0, 0.0000000000D0,    
     $       0.11142137093D0, 4.71279419371D0/     
      data stheta/ 0.88998808748D0,-0.88998808748D0, 
     $       0.39777715593D0,-0.39777715593D0, 
     $       0.86766622025D0,-0.86766622025D0 /   
      data ctheta/ 0.45598377618D0, 0.45598377618D0, 
     $       0.91748206207D0, 0.91748206207D0, 
     $       0.49714719172D0, 0.49714719172D0  /  
      data phi/  4.9368292465D0,  0.57477043300D0, 
     $       0.0000000000D0, 0.00000000000D0, 
     $       4.71279419371D0, 0.11142137093D0/



      pi = 3.1415927
      twopi   =   2.0d0*pi
      fourpi  =   4.0d0*pi
      deg_to_rad = 180.0d0/pi

c   J2000 coordinate conversions are based on the following constants
c  eps = 23.4392911111d              Obliquity of the ecliptic
c  alphaG = 192.85948d               Right Ascension of Galactic North Pole
c  deltaG = 27.12825d                Declination of Galactic North Pole
c  lomega = 32.93192d                Galactic longitude of celestial equator  
c  alphaE = 180.02322d              Ecliptic longitude of Galactic North Pole
c  deltaE = 29.811438523d            Ecliptic latitude of Galactic North Pole
c  Eomega  = 6.3839743d              Galactic longitude of ecliptic equator              


c
c if npar LT 5 then begin
c        print,' '
c        print,' 1 RA-DEC ' + equinox + ' to Galactic
c        print,' 2 Galactic       to RA-DEC' + equinox
c        print,' 3 RA-DEC ' + equinox + ' to Ecliptic'
c        print,' 4 Ecliptic       to RA-DEC' + equinox
c        print,' 5 Ecliptic       to Galactic'
c        print,' 6 Galactic       to Ecliptic'
cc
c        select = 0
c        read,'Enter selection: ',select
c endif
      I  = iselect                         ! IDL offset not for fortran
      a  = ai/deg_to_rad - phi(i)
      b = bi/deg_to_rad
      sb = sin(b) 	
      cb = cos(b)
      cbsa = cb * sin(a)
      b  = -stheta(i) * cbsa + ctheta(i) * sb
      btest = 1.0d0
      if(b.lt.btest)btest = b
      bo    = asin(btest)*deg_to_rad
c
      a =  atan2( ctheta(i) * cbsa + stheta(i) * sb, cb * cos(a))
      ao = ( mod((a+psi(i)+fourpi),twopi)) * deg_to_rad


c if ( npar EQ 2 ) then begin
c	ai = ao & bi=bo
c endif

      return
      end









        subroutine xyz(date,x,y,z)
c+
c NAME
c       XYZ
c PURPOSE:
c       Calculate geocentric X,Y, and Z  and velocity coordinates of the Sun
c EXPLANATION:
c       Calculates geocentric X,Y, and Z vectors and velocity coordinates 
c       (dx, dy and dz) of the Sun.   (The positive X axis is directed towards 
c       the equinox, the y-axis, towards the point on the equator at right 
c       ascension 6h, and the z axis toward the north pole of the equator).
c       Typical position accuracy is <1e-4 AU (15000 km).
c
c CALLING SEQUENCE:
c       XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX = ]
c
c INPUT:
c       date: reduced julian date (=JD - 2400000), scalar or vector
c
c OUTPUT:
c       x,y,z: scalars or vectors giving heliocentric rectangular coordinates
c                 (in A.U) for each date supplied.    Note that sqrt(x^2 + y^2
c                 + z^2) gives the Earth-Sun distance for the given date.
c       xvel, yvel, zvel: velocity vectors corresponding to X, Y and Z.
c
c OPTIONAL KEYWORD INPUT:
c       EQUINOX: equinox of output. Default is 1950.
c
c EXAMPLE:
c       What were the rectangular coordinates and velocities of the Sun on 
c       Jan 22, 1999 0h UT (= JD 2451200.5) in J2000 coords? NOTE:
c       Astronomical Almanac (AA) is in TDT, so add 64 seconds to 
c       UT to convert.
c
c       IDL> xyz,51200.5+64.d/86400.d,x,y,z,xv,yv,zv,equinox = 2000
c
c       Compare to Astronomical Almanac (1999 page C20)
c                   X  (AU)        Y  (AU)     Z (AU)
c       XYZ:      0.51456871   -0.76963263  -0.33376880
c       AA:       0.51453130   -0.7697110   -0.3337152
c       abs(err): 0.00003739    0.00007839   0.00005360
c       abs(err)
c           (km):   5609          11759         8040 
c
c       NOTE: Velocities in AA are for Earth/Moon barycenter
c             (a very minor offset) see AA 1999 page E3
c                  X VEL (AU/DAY) YVEL (AU/DAY)   Z VEL (AU/DAY)
c       XYZ:      -0.014947268   -0.0083148382    -0.0036068577
c       AA:       -0.01494574    -0.00831185      -0.00360365
c       abs(err):  0.000001583    0.0000029886     0.0000032077
c       abs(err)
c        (km/sec): 0.00265        0.00519          0.00557
c
c PROCEDURE CALLS:
c       PRECESS_XYZ
c REVISION HISTORY
c       Original algorithm from Almanac for Computers, Doggett et al. USNO 1978
c       Adapted from the book Astronomical Photometry by A. Henden
c       Written  W. Landsman   STX       June 1989
c       Correct error in X coefficient   W. Landsman HSTX  January 1995
c       Added velocities, more terms to positions and EQUINOX keyword,
c          some minor adjustments to calculations 
c          P. Plait/ACC March 24, 1999
c       Converted to FORTRAN!, M. Wolfire UMD for use in SIRTF GLIMPSE 
C         pipeline  February 2001
C
c-
  
c   if (n_params() eq 0) then begin
c      print,'Syntax - XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX= ]'
c      print,'     (date is REDUCED Julian date (JD - 2400000.0) )'
c      return
c   endif
      real*8 date,t,j,n,m
      pi = 3.1415927
      picon = pi/180.0d0
      t = (date - 15020.0d0)/36525.0d0   !Relative Julian century from 1900

c NOTE: longitude arguments below are given in *equinox* of date.
c   Precess these to equinox 1950 to give everything an even footing.
c   Compute argument of precession from equinox of date back to 1950
      pp = (1.396041d0 + 0.000308d0*(t + 0.5d0))*(t-0.499998d0)

c Compute mean solar longitude, precessed back to 1950
      el = 279.696678D0 + 36000.76892D0*t + 0.000303d0*t*t - pp

c Compute Mean longitude of the Moon
      c = 270.434164d0 + 480960.d0*t + 307.883142d0*t - 
     $     0.001133d0*t*t - pp

c Compute longitude of Moon's ascending node
      n = 259.183275d0 - 1800.d0*t - 134.142008d0*t + 
     $     0.002078d0*t*t - pp

c Compute mean solar anomaly
      g = 358.475833d0 + 35999.04975d0*t - 0.00015d0*t*t

c Compute the mean jupiter anomaly
      j = 225.444651d0 + 2880.0d0*t + 154.906654d0*t*t

c Compute mean anamoly of Venus
      v = 212.603219d0 + 58320.d0*t + 197.803875d0*t + 0.001286d0*t*t

c Compute mean anamoly of Mars
      m = 319.529425d0 + 19080.d0*t + 59.8585d0*t + 0.000181d0*t*t

c Convert degrees to radians for trig functions
      el = el*picon
      g  = g*picon
      j =  j*picon
      c  = c*picon
      v  = v*picon
      n  = n*picon
      m  = m*picon

c Calculate X,Y,Z using trigonometric series
      X =   0.999860d0*cos(el)                         
     $  - 0.025127d0*cos(g - el)                      
     $  + 0.008374d0*cos(g + el)                      
     $  + 0.000105d0*cos(g + g + el)                  
     $  + 0.000063d0*t*cos(g - el)                    
     $  + 0.000035d0*cos(g + g - el)                  
     $  - 0.000026d0*sin(g - el - j)                  
     $  - 0.000021d0*t*cos(g + el)                    
     $  + 0.000018d0*sin(2.d0*g + el - 2.d0*v)          
     $  + 0.000017d0*cos(c)                           
     $  - 0.000014d0*cos(c - 2.d0*el)                  
     $  + 0.000012d0*cos(4.d0*g + el - 8.d0*m + 3.d0*j)  
     $  - 0.000012d0*cos(4.d0*g - el - 8.d0*m + 3.d0*j)  
     $  - 0.000012d0*cos(g + el - v)                  
     $  + 0.000011d0*cos(2.d0*g + el - 2.d0*v)          
     $  + 0.000011d0*cos(2.d0*g - el - 2.d0*j)         
  

      Y =   0.917308d0*sin(el)                             
     $  + 0.023053d0*sin(g - el)                         
     $  + 0.007683d0*sin(g + el)                         
     $  + 0.000097d0*sin(g + g + el)                     
     $  - 0.000057d0*t*sin(g - el)                       
     $  - 0.000032d0*sin(g + g - el)                     
     $  - 0.000024d0*cos(g - el - j)                     
     $  - 0.000019d0*t*sin(g + el)                       
     $  - 0.000017d0*cos(2.d0*g + el - 2.d0*v)           
     $  + 0.000016d0*sin(c)                              
     $  + 0.000013d0*sin(c - 2.d0*el )                   
     $  + 0.000011d0*sin(4.d0*g + el - 8.d0*m + 3.d0*j)  
     $  + 0.000011d0*sin(4.d0*g - el - 8.d0*m + 3.d0*j)  
     $  - 0.000011d0*sin(g + el - v)                     
     $  + 0.000010d0*sin(2.d0*g + el - 2.d0*v )          
     $  - 0.000010d0*sin(2.d0*g - el - 2.d0*j )         


      Z =   0.397825d0*sin(el)        
     $  + 0.009998d0*sin(g-el)      
     $  + 0.003332d0*sin(g+el)     
     $  + 0.000042d0*sin(g+g+el)    
     $  - 0.000025d0*t*sin(g-el)    
     $  - 0.000014d0*sin(g+g-el)    
     $  - 0.000010d0*cos(g-el-j)    

cPrecess_to new equator?
c   if keyword_set(equinox) then precess_xyz, x, y, z, 1950, equinox
c
c   if N_params() LE 3 then return
c   
c   XVEL = -0.017200d * sin(el)           $
c          -0.000288d * sin(g + el)       $
c          -0.000005d * sin(2.d0*g + el)  $
c          -0.000004d * sin(c)            $
c          +0.000003d * sin(c - 2.d0*el)  $
c          +0.000001d *t * sin(g+el)      $
c          -0.000001d * sin(2.d0*g-el)           
c 
c   YVEL =  0.015780 * cos(el)            $
c          +0.000264 * cos(g + el)        $
c          +0.000005 * cos(2.d0*g + el)   $
c          +0.000004 * cos(c)             $
c          +0.000003 * cos(c - 2.d0*el)   $
c          -0.000001 * t * cos(g + el)    
c
c   ZVEL = 0.006843 * cos(el)             $
c         +0.000115 * cos(g  + el)        $
c         +0.000002 * cos(2.d0*g + el)    $
c         +0.000002 * cos(c)              $
c         +0.000001 * cos(c - 2.d0*el)    

cPrecess to new equator?

c   if keyword_set(equinox) then precess_xyz, xvel, yvel, zvel, 1950, equinox

      return
      end





C	SUBROUTINE INTERP
C 	
C	PURPOSE
C          INTERPOLATE BETWEEN DATA POINTS TO EVALUATE A FUNCTION
C	
C	USAGE
C	  CALL INTERP (X, Y, NPTS, NTERMS, XIN, YOUT)
C
C       DESCRIPTION OF PARAMETERS
C	  X      - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C         Y	 - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C	  NPTS   - NUMBER OF PAIRS OF DATA POINTS
C         NTERMS - NUMBER OF TERMS IN FITTING POLYNOMIAL
C         XIN    - INPUT VALUE OF X
C	  YOUT   - INTERPOLATED VALUE OF Y
C
C	SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	  NONE
C
C	MODIFICATIONS FOR FORTRAN II
C 	  OMIT DOUBLE PRECISION SPECIFICATIONS
C
C	COMMENTS
C	  DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C	  VALUE OF NTERMS MAY BE MODIFIED BY THE PROGRAM
C        * X ARRAY MUST BE IN ORDER OF INCREASING X
C
	SUBROUTINE INTERP (X, Y, NPTS, NTERMS, XIN, YOUT)
        implicit double precision (a-h,o-z)
	DOUBLE PRECISION DELTAX, DELTA, A, PROD, SUM
        DOUBLE PRECISION X,Y
	DIMENSION X(1), Y(1)
	DIMENSION DELTA(10), A(10)

C
C	  SEARCH FOR APPROPRIATE VALUE OF X(1)
C
  11  	DO 19 I=1, NPTS         
	IF (XIN-X(I)) 13, 17, 19
  13	I1 = I - NTERMS/2
	IF (I1) 15, 15, 21
  15	I1 = 1
	GO TO 21
  17	YOUT = Y(I)
  18	GO TO 61
  19	CONTINUE
	I1 = NPTS - NTERMS + 1
  21	I2 = I1 + NTERMS -1
	IF (NPTS - I2) 23, 31, 31
  23	I2 = NPTS
	I1 = I2 - NTERMS +1
  25	IF (I1) 26, 26, 31
  26 	I1 = 1
  27	NTERMS = I2 - I1 +1
C
C	  EVALUATE DEVIATIONS DELTA
C
  31	DENOM = X(I1+1) - X(I1)
	DELTAX = (XIN - X(I1)) / DENOM
	DO 35 I=1, NTERMS
	IX= I1 + I - 1
  35	DELTA(I) = (X(IX) - X(I1)) / DENOM
C
C	  ACCUMULATE COEFFICIENTS A
C
  40	A(1) = Y(I1)
  41	DO 50 K=2, NTERMS
	PROD = 1.
	SUM = 0.
	IMAX = K - 1
	IXMAX = I1 + IMAX
	DO 49 I=1, IMAX
	J = K - I
	PROD = PROD * (DELTA(K) - DELTA(J))
  49	SUM = SUM - A(J)/PROD
  50	A(K) = SUM + Y(IXMAX)/PROD
C
C	  ACCCUMULATE SUM OF EXPANSION
C
  51	SUM = A(1)
	DO 57 J=2, NTERMS
	PROD = 1.
	IMAX = J - 1
	DO 56 I=1, IMAX
  56	PROD = PROD * (DELTAX - DELTA(I))
  57	SUM = SUM + A(J)*PROD
  60	YOUT = SNGL(SUM)
  61	RETURN
	END	




