Program RAdec_Ecliptic
	implicit none
c
c	This program converts eq maps to sun-centered (or anti-sun-centered) ecliptic maps
c	This is the main program for working out the zodical light empirical model
c	The Longitude or Elongation plots here are "flipped", i.e., as seen on the sky.
c	Epoch is noon on 1 January 2000
c	Any increase in number of maps past the present 29 exceeds the memory limits of "abpc" (2 Gbytes).
c
	real*4 sind, cosd, tand, atan2d

	character*46 infile(3,29),dum
	character*5 dum2
	character*44 outfile,outfile1,outfile2	!"diff"
	character*45 outfilen
c	character*48 outfile,outfile1,outfile2	!"centered"
	character*10 head(5)
	real*4 heliomap(3600,1800,29)
	integer*4 i,j,iyr,ii,jj,isum,n,nn,mmap,ncam,mdate,nmin,nmax,k,kmax,karg,indx(29),iisv,jjsv,iii,islab,inumb
	integer*4 idoy,hh,mm,ss
	integer*1 nheliomap(3600,1800)
	integer*1 nheliomap1(3600,1800),nheliomap2(3600,1800)
	integer*4 idoy1(29),ihh1(29),imm1(29),iss1(29),idoy2(29),ihh2(29),imm2(29),iss2(29),idoy3(29),ihh3(29),imm3(29),iss3(29)
	integer*4 ixstars,iisave
	real*4 eqmap (3600,1200),doy,tconst,temp,skymap(3600,1200),skymapz(3600,1800),A8,A80,A8exp,heliomapout(3600,1800)!,cPA,cPA8
	real*4 rLng,rLngi,rLat,rDis,rLngSun,rrLng,rrLat,sum,array(29),gain,c,s,zod,ssq,arg,aEloP,rrrLng,rrrLat,zslab,argLat,fractNS
	real*8 dLngSun,dLatSun,dDisSun
	real*4 theta1,rEloP,rP,PA,R,EloSun,Omega,adhocNS,adhocEW,RRR,REW,RUSE,dis(29)
	real*4 heliomap1(3600,1800),heliomap2(3600,1800)
	real*4 arg1,arg2,arg3,argx,argy,weight
	data dum2 /'Orbit'/
c
	data ixstars /1/ !ixstars=1 means: yes, use removed-star files in EngineeringXStars
c	
	data tconst /6096.312/
c
c	data R,Omega,RRR,REW /3.,90.,3.,12./			!values 6  Nov 2006	!1st 2 values from AQ, 4th, page 330
c	data R,Omega,RRR,REW /5.,70.,6.,6./			!values 16 Nov 2006	!Some evidence for a smaller Omega...
c	data R,Omega,RRR,REW /5.,90.,6.,0./			!values October 2008, ***REW=0 kills EW asymmetry everywhere***
	data R,Omega,RRR,REW /5.,78.25,6.,0./			!values 14 November 2008: note 0.75*value RRR is used in North, CBZodi value for Omega
	data A8,A80,A8exp / 8800.,1200.,10./			!values 18 November 2008: PA dependent term; an old name, but now exponential in 2D	
	data gain /1.0/
c	
	data idoy1 /56,127,183,242,319,361,55,113,172,238,291,345,41,101,184,247,304,364,61,149,204,264,323, 76,138,194, 06,125,184/
	data ihh1  /16, 15, 11, 15, 15, 01,12, 16, 17, 16, 01, 00,00, 00, 01, 01, 01, 00,00, 00, 00, 00, 01, 23, 03, 00, 00, 01, 05/
	data imm1  /05, 41, 10, 56, 25, 25,54, 12, 30, 44, 10, 31,54, 10, 06, 10, 17, 30,50, 14, 55, 07, 18, 57, 55, 53, 56, 13, 44/
	data iss1  /14, 11, 38, 25, 03, 01,03, 13, 21, 44, 02, 41,32, 30, 55, 04, 09, 52,08, 43, 32, 32, 19, 58, 30, 23, 51, 37, 23/
c	
	data idoy2 /58,124,182,240,311,358,53,112,171,241,290,346,38,100,185,248,305, 01,56,148,206,266,325, 79,137,196, 07,127,186/
	data ihh2  /15, 16, 04, 16, 16, 15,01, 16, 17, 04, 01, 00,00, 00, 00, 00, 00, 01,00, 00, 00, 01, 00, 01, 05, 00, 00, 00, 05/
	data imm2  /30, 33, 41, 31, 04, 50,38, 29, 47, 00, 27, 14,05, 28, 49, 52, 59, 37,37, 32, 20, 13, 42, 04, 54, 17, 39, 38, 08/
	data iss2  /14, 46, 45, 31, 01, 35,03, 50, 58, 40, 41, 02,56, 11, 14, 22, 27, 01,07, 27, 05, 39, 50, 03, 50, 51, 04, 02, 48/
c	
	data idoy3 /60,129,194,198,238,316,356,50,111,170,240,296,347,39, 97,187,246,306,365,58,147,205,265,324,078,136,195,325, 08/
	data ihh3  /16, 16, 01, 00, 15, 16, 16,17, 15, 18, 04, 14, 01,01, 06, 00, 01, 04, 00,00, 00, 00, 01, 01, 01, 04, 00, 17, 00/
	data imm3  /36, 47, 11, 01, 25, 17, 25,45, 05, 05, 18, 56, 37,29, 25, 13, 27, 04, 13,01, 50, 37, 31, 00, 21, 31, 35, 57, 21/
	data iss3  /51, 44, 18, 09, 01, 46, 46,15, 51, 35, 18, 13, 59,52, 58, 51, 45, 56, 09,41, 10, 48, 23, 35, 48, 01, 37, 47, 18/
c
	data kmax /29/	
c	data outfile     /'f:\EngineeringXStars\caeq2_centered_longlata.grd'/
c	data outfile1    /'f:\EngineeringXstars\caeq2_centered_pa__epsa.grd'/
c	data outfile2    /'f:\EngineeringXstars\caeq2_centered_eps_lata.grd'/
	data outfile     /'f:\EngineeringXStars\caeq2_diff_longlata.grd'/
	data outfilen    /'f:\EngineeringXStars\caeq2_diff_longlatBn.grd'/
c	data outfile1    /'f:\EngineeringXstars\caeq2_diff_pa__epsa.grd'/
	data outfile2    /'f:\EngineeringXstars\caeq2_diff_eps_lata.grd'/
c	
	data infile(1,1) /'f:\EngineeringXStars\c1eq2_2003_056_160514.grd'/
	data infile(1,2) /'f:\EngineeringXStars\c1eq2_2003_127_154111.grd'/
	data infile(1,3) /'f:\EngineeringXStars\c1eq2_2003_183_111038.grd'/
	data infile(1,4) /'f:\EngineeringXStars\c1eq2_2003_242_155625.grd'/
	data infile(1,5) /'f:\EngineeringXStars\c1eq2_2003_319_152503.grd'/
	data infile(1,6) /'f:\EngineeringXStars\c1eq2_2003_361_012501.grd'/
	data infile(1,7) /'f:\EngineeringXStars\c1eq2_2004_055_125403.grd'/
	data infile(1,8) /'f:\EngineeringXStars\c1eq2_2004_113_161213.grd'/
	data infile(1,9) /'f:\EngineeringXStars\c1eq2_2004_172_173021.grd'/
	data infile(1,10) /'f:\EngineeringXStars\c1eq2_2004_238_164444.grd'/
	data infile(1,11) /'f:\EngineeringXStars\c1eq2_2004_291_011002.grd'/
	data infile(1,12) /'f:\EngineeringXStars\c1eq2_2004_345_003141.grd'/
	data infile(1,13) /'f:\EngineeringXStars\c1eq2_2005_041_005432.grd'/
	data infile(1,14) /'f:\EngineeringXStars\c1eq2_2005_101_001030.grd'/
	data infile(1,15) /'f:\EngineeringXStars\c1eq2_2005_184_010655.grd'/
	data infile(1,16) /'f:\EngineeringXStars\c1eq2_2005_247_011004.grd'/
	data infile(1,17) /'f:\EngineeringXStars\c1eq2_2005_304_011709.grd'/
	data infile(1,18) /'f:\EngineeringXStars\c1eq2_2005_364_003052.grd'/
	data infile(1,19) /'f:\EngineeringXStars\c1eq2_2006_061_005008.grd'/
	data infile(1,20) /'f:\EngineeringXStars\c1eq2_2006_149_001443.grd'/
	data infile(1,21) /'f:\EngineeringXStars\c1eq2_2006_204_005532.grd'/
	data infile(1,22) /'f:\EngineeringXStars\c1eq2_2006_264_000732.grd'/
	data infile(1,23) /'f:\EngineeringXStars\c1eq2_2006_323_011819.grd'/
	data infile(1,24) /'f:\EngineeringXStars\c1eq2_2007_076_235758.grd'/
	data infile(1,25) /'f:\EngineeringXStars\c1eq2_2007_138_035530.grd'/
	data infile(1,26) /'f:\EngineeringXStars\c1eq2_2007_194_005323.grd'/
	data infile(1,27) /'f:\EngineeringXStars\c1eq2_2008_006_005651.grd'/
	data infile(1,28) /'f:\EngineeringXStars\c1eq2_2008_125_011337.grd'/
	data infile(1,29) /'f:\EngineeringXStars\c1eq2_2008_184_054423.grd'/
c	
	data infile(2,1) /'f:\EngineeringXStars\c2eq2_2003_058_153014.grd'/
	data infile(2,2) /'f:\EngineeringXStars\c2eq2_2003_124_163346.grd'/
	data infile(2,3) /'f:\EngineeringXStars\c2eq2_2003_182_044145.grd'/
	data infile(2,4) /'f:\EngineeringXStars\c2eq2_2003_240_163131.grd'/
	data infile(2,5) /'f:\EngineeringXStars\c2eq2_2003_311_160401.grd'/
	data infile(2,6) /'f:\EngineeringXStars\c2eq2_2003_358_155035.grd'/
	data infile(2,7) /'f:\EngineeringXStars\c2eq2_2004_053_013803.grd'/
	data infile(2,8) /'f:\EngineeringXStars\c2eq2_2004_112_162950.grd'/
	data infile(2,9) /'f:\EngineeringXStars\c2eq2_2004_171_174758.grd'/
	data infile(2,10) /'f:\EngineeringXStars\c2eq2_2004_241_040040.grd'/
	data infile(2,11) /'f:\EngineeringXStars\c2eq2_2004_290_012741.grd'/
	data infile(2,12) /'f:\EngineeringXStars\c2eq2_2004_346_001402.grd'/
	data infile(2,13) /'f:\EngineeringXStars\c2eq2_2005_038_000556.grd'/
	data infile(2,14) /'f:\EngineeringXStars\c2eq2_2005_100_002811.grd'/
	data infile(2,15) /'f:\EngineeringXStars\c2eq2_2005_185_004914.grd'/
	data infile(2,16) /'f:\EngineeringXStars\c2eq2_2005_248_005222.grd'/
	data infile(2,17) /'f:\EngineeringXStars\c2eq2_2005_305_005927.grd'/
	data infile(2,18) /'f:\EngineeringXStars\c2eq2_2006_001_013701.grd'/
	data infile(2,19) /'f:\EngineeringXStars\c2eq2_2006_056_003707.grd'/
	data infile(2,20) /'f:\EngineeringXStars\c2eq2_2006_148_003227.grd'/
	data infile(2,21) /'f:\EngineeringXStars\c2eq2_2006_206_002005.grd'/
	data infile(2,22) /'f:\EngineeringXStars\c2eq2_2006_266_011339.grd'/
	data infile(2,23) /'f:\EngineeringXStars\c2eq2_2006_325_004250.grd'/
	data infile(2,24) /'f:\EngineeringXStars\c2eq2_2007_079_010403.grd'/
	data infile(2,25) /'f:\EngineeringXStars\c2eq2_2007_137_055450.grd'/
	data infile(2,26) /'f:\EngineeringXStars\c2eq2_2007_196_001751.grd'/
	data infile(2,27) /'f:\EngineeringXStars\c2eq2_2008_007_003904.grd'/
	data infile(2,28) /'f:\EngineeringXStars\c2eq2_2008_127_003802.grd'/
	data infile(2,29) /'f:\EngineeringXStars\c2eq2_2008_186_050848.grd'/
c
	data infile(3,1) /'f:\EngineeringXStars\c3eq2_2003_060_163651.grd'/
	data infile(3,2) /'f:\EngineeringXStars\c3eq2_2003_129_164744.grd'/
	data infile(3,3) /'f:\EngineeringXStars\c3eq2_2003_194_011118.grd'/
	data infile(3,4) /'f:\EngineeringXStars\c3eq2_2003_198_000109.grd'/
	data infile(3,5) /'f:\EngineeringXStars\c3eq2_2003_238_152501.grd'/
	data infile(3,6) /'f:\EngineeringXStars\c3eq2_2003_316_161746.grd'/
	data infile(3,7) /'f:\EngineeringXStars\c3eq2_2003_356_162546.grd'/
	data infile(3,8) /'f:\EngineeringXStars\c3eq2_2004_050_174515.grd'/
	data infile(3,9) /'f:\EngineeringXStars\c3eq2_2004_111_150551.grd'/
	data infile(3,10) /'f:\EngineeringXStars\c3eq2_2004_170_180535.grd'/
	data infile(3,11) /'f:\EngineeringXStars\c3eq2_2004_240_041818.grd'/
	data infile(3,12) /'f:\EngineeringXStars\c3eq2_2004_296_145613.grd'/
	data infile(3,13) /'f:\EngineeringXStars\c3eq2_2004_347_013759.grd'/
	data infile(3,14) /'f:\EngineeringXStars\c3eq2_2005_039_012952.grd'/
	data infile(3,15) /'f:\EngineeringXStars\c3eq2_2005_097_062558.grd'/
	data infile(3,16) /'f:\EngineeringXStars\c3eq2_2005_187_001351.grd'/
	data infile(3,17) /'f:\EngineeringXStars\c3eq2_2005_246_012745.grd'/
	data infile(3,18) /'f:\EngineeringXStars\c3eq2_2005_306_040456.grd'/
	data infile(3,19) /'f:\EngineeringXStars\c3eq2_2005_365_001309.grd'/
	data infile(3,20) /'f:\EngineeringXStars\c3eq2_2006_058_000141.grd'/
	data infile(3,21) /'f:\EngineeringXStars\c3eq2_2006_147_005010.grd'/
	data infile(3,22) /'f:\EngineeringXStars\c3eq2_2006_205_003748.grd'/
	data infile(3,23) /'f:\EngineeringXStars\c3eq2_2006_265_013123.grd'/
	data infile(3,24) /'f:\EngineeringXStars\c3eq2_2006_324_010035.grd'/
	data infile(3,25) /'f:\EngineeringXStars\c3eq2_2007_078_012148.grd'/
	data infile(3,26) /'f:\EngineeringXStars\c3eq2_2007_136_043101.grd'/
	data infile(3,27) /'f:\EngineeringXStars\c3eq2_2007_195_003537.grd'/
	data infile(3,28) /'f:\EngineeringXStars\c3eq2_2007_325_175747.grd'/
	data infile(3,29) /'f:\EngineeringXStars\c3eq2_2008_008_002118.grd'/
c
	data heliomap /187920000*0./	!29 maps
	data heliomap1 /6480000*0./
	data heliomap2 /6480000*0./
	data nheliomap /6480000*0/
	data nheliomap1 /6480000*0/
	data nheliomap2 /6480000*0/
	data heliomapout /6480000*0./
	data mmap /2/	!mmap=0 means skip outputing map
			!mmap=1 means output only Long-Lat map
			!mmap=2 means output only PA-Eps   map
			!mmap=3 means output only Eps-Lat  map
			!mmap=4 means output all 3 maps
c
c	mdate=0: use fixed date (Jan 1.5, 2000)
c	mdate=1: Place Sun at 180 degrees in Long-Lat Plot
c	mdate=2: Sun split at 0, 360 degrees in Long-Lat Plot
c	
	data mdate /1/	
	data nmin,nmax / 1,29 /
	data outfile1     /'f:\EngineeringXStars\caeq2_dif_p_epc3_04.grd'/
c
	data islab / 1 /	!islab=0 means skip adding ad hoc slab
	data inumb /0/		!inumb=1 means yes, output the number distribution
c
c	open(30,file='f:\EngineeringXStars\tempdist.dat')	!This file was used to output antisolar square results
c		
	if(ixstars.eq.0)then
		do nn=1,3
		do n=1,29
			dum=infile(nn,n)
			write(dum(15:19),'(a5)')dum2
			infile(nn,n)=dum
		enddo
		enddo
		write(outfile(15:19),'(a5)')dum2
		write(outfile1(15:19),'(a5)')dum2
		write(outfile2(15:19),'(a5)')dum2
	endif
c
c	Get skymaps to be subtracted
c
	print *,'reading skymaps'
	if(ixstars.eq.0)open(10,file='f:\EngineeringOrbits\eq_allskymeed2.grd',status='old')
c	if(ixstars.eq.1)open(10,file='f:\EngineeringXStars\eq_allskymeed2.grd',status='old')
	if(ixstars.eq.1)open(10,file='f:\EngineeringXStars\eq_allskymeed_avg12.grd',status='old')
	read(10,11)head
	do j=1,1200
		read (10,110)(skymap(i,j),i=1,3600)
	enddo
	close(10)
c
c	This map, used below to discard way-off data 2X too large, should be updated when model development is complete, won't change much.
c	The one currently in use comes from 10/3/2006, seems to work pretty well. A newer one is more noisy, don't know why...
c	It does a pretty good job of reducing noise around subtracted bright stars or planets, and bright shutter edges for camera 3.
c	This map is made with the 10.*log10 option below, without zod subtracted or in ratio, then just change the name, as indicated below.
c	Strictly, this map ought to be remade each time something is changed in "parameter hell", but the cut is so crude it probably matters but little.
c
	if(ixstars.eq.0)open(11,file='f:\EngineeringOrbits\caeq2_centered_longlat_all.grd',status='old')!then rename the output map with "_all" to get final
	if(ixstars.eq.1)open(11,file='f:\EngineeringXStars\caeq2_centered_longlat_all.grd',status='old')!then rename the output map with "_all" to get final
	read(11,11)head
	do j=1,1800
		read (11,110)(skymapz(i,j),i=1,3600)
	enddo
	close(11)
c
c	do nn=3,3	!Only camera 1
c	do nn=2,2	!Only camera 2
	do nn=1,1	!Only camera 3
c	do nn=1,3	!The next statement selects the order in which one camera "rules" over the next, later camera rules when more than one...
	ncam=4-nn	!This establishes the heirarchy cameras 1, then 2, then 3, which is the one Bernie uses.
c	ncam=nn		!This establishes the heirarchy cameras 3, then 2, then 1, opposite.
	print *,'starting camera ',ncam
	do n=nmin,nmax
c****************		time-of-year runs, cameras 1&2
c		if(n.ne.1.and.n.ne.7.and.n.ne.13.and.n.ne.19.and.n.ne.24)go to 95	!1st period
c		if(n.ne.2.and.n.ne.8.and.n.ne.14.and.n.ne.28.and.n.ne.25)go to 95	!2nd period
c		if(n.ne.3.and.n.ne.9.and.n.ne.15.and.n.ne.20.and.n.ne.29)go to 95	!3rd period
		if(n.ne.4.and.n.ne.10.and.n.ne.16.and.n.ne.21.and.n.ne.26)go to 95	!4th period
c		if(n.ne.5.and.n.ne.11.and.n.ne.17.and.n.ne.22)go to 95			!5th period
c		if(n.ne.6.and.n.ne.12.and.n.ne.18.and.n.ne.23.and.n.ne.27)go to 95	!6th period
c****************	
c		***Note camera 3, n=18 (2005,DOY306) has bad quaternions near G.C. Thought to be due to nearby Venus
c
		if(mdate.eq.0)then
c
c        	  JD2000 = 2451545.D0       ! Jan 1.5, 2000 (noon), put this date in for all maps, if mdate=0
c
	  	  iyr=2000
		  doy=1.5
		else	!calculate appropriate date
c
c		  Insert year
c
		  iyr=2003
	if(ncam.eq.1)then
		if(n.gt.6)iyr=iyr+1
		if(n.gt.12)iyr=iyr+1
		if(n.gt.18)iyr=iyr+1
		if(n.gt.23)iyr=iyr+1
		if(n.gt.26)iyr=iyr+1
	endif
	if(ncam.eq.2)then
		if(n.gt.6)iyr=iyr+1
		if(n.gt.12)iyr=iyr+1
		if(n.gt.17)iyr=iyr+1
		if(n.gt.23)iyr=iyr+1
		if(n.gt.26)iyr=iyr+1
	endif
	if(ncam.eq.3)then
		if(n.gt.7)iyr=iyr+1
		if(n.gt.13)iyr=iyr+1
		if(n.gt.19)iyr=iyr+1
		if(n.gt.24)iyr=iyr+1
		if(n.gt.28)iyr=iyr+1
	endif
c
c		  Insert idoy, hour, minutes, seconds		
c
		  if(ncam.eq.1)then
			idoy=idoy1(n)
			hh=ihh1(n)
			mm=imm1(n)
			ss=iss1(n)
		  endif
		  if(ncam.eq.2)then
			idoy=idoy2(n)
			hh=ihh2(n)
			mm=imm2(n)
			ss=iss2(n)
		  endif
		  if(ncam.eq.3)then
			idoy=idoy3(n)
			hh=ihh3(n)
			mm=imm3(n)
			ss=iss3(n)
		  endif
		  doy=float(idoy)+(3600.*float(hh)+60.*float(mm)+float(ss)+0.5*tconst)/86400.
c	
c		  shift zodiac calculation to middle of time window, sort of...	
c
		  if(iyr.eq.2003.or.(iyr.eq.2004.and.idoy.lt.100).or.(iyr.eq.2006.and.idoy.eq.56))then
		    doy=doy+1.
		  else
		    doy=doy+0.5
		  endif
c
c		Add a bit more for camera 2 2nd file, take away for camera 2 3rd file & camera 1 16th file, and add for camera 3 4th file.
c
		if(ncam.eq.2.and.n.eq.2)doy=doy+0.5
		if(ncam.eq.2.and.n.eq.3)doy=doy-0.35
		if(ncam.eq.1.and.n.eq.16)doy=doy-0.35
		if(ncam.eq.3.and.n.eq.4)doy=doy+0.2
c
		endif	!endif for the fixed date throughout, haven't used this for a long time, forgotten why it's included as an option here...
c
c		The fitting of stars requires 1.6%/yr to flatten the response, but looking at Galactic Center requires camera 2 @ 1.0%/yr.
c		A look at camera 1 Galactic Center requires 0.5%/yr.
c		
		if(ncam.gt.1)gain=1.00+0.01*(float(iyr-2003)+doy/365.25) !gain correction for CCD response loss vs time: note we use 1%/yr here, not 1.6%/yr.
		if(ncam.eq.1)gain=1.00+0.01*(float(iyr-2003)+doy/365.25) !gain correction for CCD response loss vs time: note we use 1%/yr here, not 1.6%/yr.
		if(ncam.eq.1)gain=gain*0.97			!camera relative gains, Chris Eyles feels these #s should be much closer to unity
		if(ncam.eq.3)gain=gain*0.92			!Camera 3 Value for 16 Feb 06 params...
c		
		call SunNewcomb(0,iyr,doy,dLngSun,dLatSun,dDisSun)
		rDis=sngl(dDisSun)
		dis(n)=rDis
		rLng=sngl(dLngSun)
		rLngSun = rLng
		rrLng=rLng
		rLat=sngl(dLatSun)
		rrLat=rLat
		adhocNS=cosd(rLngSun-Omega-90.)
		adhocEW=cosd(rLngSun-Omega)
		call ECLIPTIC_EQUATOR(0,iyr,doy,rrLng,rrLat)
		print *,n,'Sun at Ecliptic Long',rLng,', DOY',doy,': RA,Dec =',rrLng,rrLat
c		print *,n,'Sun-Earth distance = ',rDis
		print *,n,'adhocNS = ',adhocNS
c		
		open(11,file=infile(ncam,n),status='old')
		read (11,11)head
11		format(a10)
		do j=1,1200
			read (11,110)(eqmap(i,j),i=1,3600)
110			format(20f8.1)	
		enddo
		close(11)
c
	do i=1,3600
      	    rLngi=360.-(float(i)/10.-0.05)	!correct for flipped map
	    do j=1,1200
	    	rLng=rLngi
		rLat=-60.05+float(j)/10.
		call ECLIPTIC_EQUATOR(1,iyr,doy,rLng,rLat)
		if(mdate.eq.0)ii=nint(10.*(rLng+0.05))
		if(mdate.eq.1)ii=nint(10.*(rLng-rLngSun-180.+0.05))
		if(mdate.eq.2)ii=nint(10.*(rLng-rLngSun+0.05))
		if(ii.lt.1)ii=ii+3600
		if(ii.lt.1)ii=ii+3600		!sometimes you need two, because of the 180 deg offset
		if(ii.gt.3600)ii=ii-3600	!shouldn't ever happen
		if(ii.lt.1)ii=1			!shouldn't ever happen
		ii=3601-ii			!flipped map restore
		iisv=ii
		iisave=ii
		if(mdate.eq.2)then
			iisave=iisave-1800	!could as well be "+1800", if sign flipped next statement...
			if(iisave.lt.0)iisave=iisave+3600
		endif
		jj=nint(10.*(rLat+90.05))
		jjsv=jj
c
c		Impose relative camera gain, subtract Sidereal Sky using "skymap".
c		Then correct for distance to Sun using AQ value = rDis**2.3.
c
		if(eqmap(i,j).gt.0.)then			!non-zero map entries continue on...
		    eqmap(i,j)=eqmap(i,j)*gain-skymap(i,j)
		    eqmap(i,j)=eqmap(i,j)*(rDis**2.3)
		    arg=10.*alog10(eqmap(i,j))
		    if(mdate.eq.0.or.arg.lt.(skymapz(iisave,jj)+3.0))then	!allow contributions up to 2X of naive "0" or cleaned "_a" average 
			if(jj.gt.1800.or.jj.lt.1)print *,'1',jj		
c			
			rLng=rLng-rLngSun			!Geocentric ecliptic
			if(rLng.le.0.)rLng=rLng+360.		!probably not needed...
			rP=max(0.25,cosd(rLng)*cosd(rLat))	!Geocentric distance
			PA=atan2d(sind(rLng),tand(rLat))
			rrrLng=rLng
			rrrLat=rLat
			call POINT_ON_LOS(rrrLng,rrrLat,rP,rEloP,EloSun,iii)	!iii=-1 is West
c@#$%!			
c			if(rEloP.lt.60.)go to 94	!@#$%!  Limit elongation range for making siderial maps...(not done here)
c@#$%!			
			if(rEloP.lt.0.)print *,'negative rEloP = ',rEloP
			theta1=PA		!+R*cosd(rLngSun-Omega)	!*****this final term does little, ignore it
c			cPA=sind(theta1)
c			cPA8=cPA**8
c			
			aEloP=180.-rEloP
			c=cosd(rEloP)
			ssq=1.-c*c
			s=sqrt(ssq)
c
c			Here set up altering the angular scale in like fashion to the above (rDis**2.3) correction, see comment just below
c			
			ssq=s**2.3
c
c	****Entering Parameter Hell**** 
c	{"Abandon all understanding, ye who enter here..."}
c
c	Calculate expected zodiacal-light contribution, using the first of Steve Price's two suggested functional forms.
c
c	This is based on Lamy and Perrin (Astr. Ap. 163, 269 (1986)), with a VSF from Ney and Merrill (Science 194, 1051 (1976)),
c	which integrates to a sum of cosines up through power 3, all divided by (rDis sin(rEloP))**2.3.
c				
c			zod=75.	+120.*c		!values 15 Nov 2006
c			zod=43.	+120.*c		!values 15 April 2008 with full gaussian below
			zod=65.	+120.*c		!values 28 April 2008 **********************
			if(c.le.0.)then		!antisolar hemisphere
				zod=zod+154.*c*c+88.*(c**3)				!values 17 Nov 2006; don't forget c**3 is (-) here...
			else			!solar hemisphere
c				zod=zod-185.*c*c+140.*(c**3)	!values 15 Nov 2006
				zod=zod-185.*c*c+65.*(c**3)	!values 18 Nov 2008, reduced to compensate for added PA term below
				zod=zod/ssq
			endif
c
c	Dropoff exponential in sind(rLat), with rEloP slope
c	The cusp is removed by using a hyperbola in place of rLat in the exponential
c	This has a flat bottom, pulls the rest of the curve down by 1.5 degrees, with the transition at about 1.5 degrees
c	This may generate another cusp at +-90  deg., hopefully too small to see...{to fix, multiply by 90/88.5}
c	The rEloP offset below comes from observed slopes on the rLat - rEloP plot of log10 (Observed/Model)
c	An offset range between about 30 to 45 is permissable.
c	
			argLat=1.5*(sqrt(1.+0.444*rLat*rLat)-1.)		!hyperbola
c			arg=abs(sind(argLat))/(0.007*(rEloP+40.))		! 2 Nov   2006 values
			arg=abs(sind(argLat))/(0.009*(rEloP+40.))		!28 April 2008 values
			zod=zod*(10**(-arg))
c
c	Ad hoc Latitude adjustments
c
			zod=zod+8.*(1.0-cosd(argLat))				!23 June 2008, see also dc offset below
c
c	For EloP < 90 degrees, include another contribution dropping off as cos(rLat).
c
			if(rEloP.lt.90.)then
				zod=zod+(1./ssq - 1.)*30.*cosd(argLat)	!16 Nov 2006
			endif
c
c	This gaussian contribution is used with camera 2 when making C2{eq,np,sp}3 grid files to be used in allskymap.for. It is
c	empirically found to clear out a broad ecliptic band in the eq siderial map, so the siderial residue has "no memory" of the ecliptic.
c	Here, "no memory" means down below a few ADUs, as best we can see...
c	This arbitrarily defines the "siderial component" to be subtracted from all maps, and "if it looks good, it's good".
c	It doesn't strictly constrain the zodiacal-model parameters, although we should probably minimize subsequent changes. 
c
			zod=zod+ 6.*exp(-rLat*rLat/512.)	!gaussian contribution sigma = 16 deg., 19   May 2008, works best for camera 1
c
c	Here subtract the double-exponential Gegenschein terms, as in our Icarus manuscript
c	The weight formula here is an empirical match to the observed Gegenschein distribution near the antisolar direction.
c	It uses 2% of the sine^2 of twice the angle times the antisolar elongation (deg.) to diminish the geometric weighting of the x,y components.
c	This is ugly, but reproduces the lozenge shape that the Gegenschein has.
c	This formula gets somewhat crazy for aEloP > 90 degrees, so we reduce Gegenschein to < 0.1 ADU at aEloP > 90 deg., using a Gaussian rolloff
c
			argx=(180.-abs(rLng))/aEloP
			argy=rLat/aEloP
			arg1=7.5*exp(-aEloP/4.)+39.5*exp(-aEloP/25.)	!This one for Latitude
			arg2=7.5*exp(-aEloP/4.)+39.5*exp(-aEloP/35.)	!This one for Longitude
			arg3=argy*argy*arg1+argx*argx*arg2
			weight=(argx*argy)**2
			weight=1.-0.02*weight*aEloP
			if(weight.lt.0.)weight=0.			!Kill negative values			
			if(aEloP.gt.60.)weight=weight*exp(-((aEloP-60.)**2)/300.)	!Roll it down fast for aEloP > 60 degrees
			zod=zod+weight*arg3
c
c	Here add a position-angle term which drops off exponentially with distance from the Sun, thus negligible in the antisolar hemisphere.
c
c			zod=zod+A8*cPA8*exp(-rEloP/A8exp)	!April 2008. Old. This varies like the 8th power of the sine of the position angle.
c
c	Here use a double exponential, with a 3 degree hyperbola derived from the deviation of PA from 90 degrees.
c
			arg=abs(PA)-90.
			zod=zod+(A8*exp(-(sqrt(1.+0.111*arg*arg)-1.)/10.)-A80)*exp(-rEloP/A8exp)	!18 November 2008
c			
c	The ad hoc dust slab...
c
			if(islab.ne.0)then
			   zslab=0.
			   fractNS=1.
			   RUSE=RRR
			   if(rEloP.gt.90.)then
				RUSE=RUSE*s		!Drop off the slab brightness in antisolar hemisphere
			   else
		  		RUSE=RUSE/ssq		!Increase slab brightness towards Sun.		
c		  		fractNS=s
			   endif
c 		        			
			   if(adhocNS.lt.0.)then	!Ad hoc N for dust slab above
			   	if(rElop.gt.90.)then
			   		RUSE=0.75*RUSE		!Use 3/4 as much in the North in antisolar hemisphere
			   	else
			   	   	RUSE=(0.75+2.*c)*RUSE	!Use 3/4 as much in the North, but more toward Sun.
				endif
				if(rLat.gt.R)zslab=-RUSE*adhocNS
				if(rLat.lt.R.and.rLat.gt.-R)zslab=-RUSE*adhocNS*(R+rLat)/(2.*R)	!blur edge by +-R
			   endif
			   if(adhocNS.gt.0.)then	!Ad hoc S for dust slab below
			   	if(rElop.le.90.)RUSE=(1.+2.0*c)*RUSE	!More towards Sun in the South.
				if(rLat.lt.-R)zslab=RUSE*adhocNS
				if(rLat.lt.R.and.rLat.gt.-R)zslab=RUSE*adhocNS*(R-rLat)/(2.*R)
			   endif
c
c	Close to the Sun the NS asymmetry is less and we include a symmetric slab
c			   
c			   zslab=zslab*fractNS+(1.-fractNS)*RUSE*abs(adhocNS)
c
c	Here's a ugly attempt at the EW/NS component, greatest when Earth passes through rotated dust plane.
c	There's got to be a better way than this: maybe building and integrating a 3-dimensional model isn't so bad...
c
c			   arg=abs(R*adhocEW*sind(rEloP))	!bend it by the R degree tilt...
c			   RUSE=REW*abs(c)			!cosine added 9 Nov to prevent E-W asymmetry, ugly...
c			   if(rEloP.gt.90.)then
c				RUSE=RUSE*s		!Drop off the slab brightness in antisolar hemisphere
c			   else
c		  		RUSE=RUSE/ssq		!Increase slab brightness towards Sun.
c			   endif
c
c	Here the slab is placed on the sky. Feathered along the equatorial line by the R line tilted.
c	Feathered exponentially across the NS line, near zero and 180 longitude, by -longitude/R/2 antisolar and by /3 solar
c			   
c			   if(adhocEW.lt.0.)then
c				if(iii.eq.-1)then	!West
c					if(rLat.ge.0.)zslab=zslab-RUSE*adhocEW
c					if(rLat.lt.0..and.rLng.gt.-180.)zslab=zslab-RUSE*adhocEW*exp(-(180.+rLng)*2./R)	!feather NS line using R/2 dropoff
c					if(rLat.lt.0..and.rLng.lt.0.)zslab=zslab-RUSE*adhocEW*exp(rLng*3./R)		!feather NS line using R/3 dropoff
c					if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab-RUSE*adhocEW*(arg+rLat)/arg		!feather equatorial band
c				endif
c				if(iii.eq.+1)then	!East
c					if(rLat.le.0.)zslab=zslab-RUSE*adhocEW
c					if(rLat.gt.0..and.rLng.lt.-180.)zslab=zslab-RUSE*adhocEW*exp((180.+rLng)*2./R)	!feather NS line using R/2 dropoff
c					if(rLat.gt.0..and.rLng.gt.0.)zslab=zslab-RUSE*adhocEW*exp(-rLng*3./R)		!feather NS line using R/3 dropoff
c					if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab-RUSE*adhocEW*(arg-rLat)/arg		!feather equatorial band
c				endif
c			   else
c				if(iii.eq.-1)then	!West
c					if(rLat.le.0.)zslab=zslab+RUSE*adhocEW
c					if(rLat.gt.0..and.rLng.gt.180.)zslab=zslab+RUSE*adhocEW*exp((180.-rLng)*2./R)	!feather NS line using R/2 dropoff
c					if(rLat.gt.0..and.rLng.lt.0.)zslab=zslab+RUSE*adhocEW*exp(rLng*3./R)		!feather NS line using R/3 dropoff
c					if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab+RUSE*adhocEW*(arg-rLat)/arg		!feather equatorial band
c				endif
c				if(iii.eq.+1)then	!East
c					if(rLat.ge.0.)zslab=zslab+RUSE*adhocEW
c					if(rLat.lt.0..and.rLng.lt.180.)zslab=zslab+RUSE*adhocEW*exp(-(180.-rLng)*2./R)	!feather NS line using R/2 dropoff
c					if(rLat.lt.0..and.rLng.gt.0.)zslab=zslab+RUSE*adhocEW*exp(-rLng*3./R)		!feather NS line using R/3 dropoff
c					if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab+RUSE*adhocEW*(arg+rLat)/arg		!feather equatorial band
c				endif
c			   endif
c
			   zod=zod+zslab
			endif
c			
c	note when rEloP=90 and rLat=90 (ecliptic poles), this + DC offset below should agree with AQ p329 (22.4), or p330 (24)
c
c	Nov 2 2006's value is 6 {from the 75 constant above top line for zod, times 10**(-1/(.007*130))}, plus the value below; but see comment just below			
c	28 April's 2008 value is 9.1 {from the 65 constant above top line for zod, times 10**(-1/(.009*130))}; but see just below
c	Then we add another 8 for the latitude adjustment, and finally the below number for the comparison, works out OK.
c
c			zod=zod+16.	!Nov 2006 constant offset, a case can be made that this should be "zod=zod+REW", but unclear how to include the slab
			if(ncam.eq.1)zod=zod+7.	!value  7 Nov. 2008
			if(ncam.eq.2)zod=zod+7.	!value 29 Oct. 2008
			if(ncam.eq.3)zod=zod-3.	!value 25 Nov. 2008
c
c	Diminish zodical light a bit near ecliptic and around 90 degrees, this probably needs refinement, especially towards the Sun.
c	This was added to take account of the "dust bands", which were created here by punching down the middle of the ecliptic -- hokey.
c			
c			zod=zod-8.*sind(rEloP)*exp(-rLat*rLat/10.)	!added 17 Nov 2006, commented out April 2008
c			
c
c	****Leaving Parameter Hell****
c
			temp=eqmap(i,j)-zod	!differences
c
c			Here put a cut which reduces moon for camera 1, 2004, DOY 238
c
			if(ncam.eq.1.and.n.eq.10.and.temp.gt.15.)go to 94
c**************
c			Note at this place in the program "zodiac_modelb.for", a distance correction is made, but only to zod:
c			temp=temp/(rDis**2.3)	!Correction to 1 AU: heliomap is reduced for > 1 AU and v.v.
c			Here, the correction is done to eqmap(i,j) above...probably doesn't matter much, but inconsistent...
c			The choice here is to show exactly the deviations SMEI sees, while elsewhere results are normalized to 1 AU.
c**************
c
c			temp=eqmap(i,j)		!direct map, use this when making the caeq2_centered_longlat_all.grd to be used above
c			temp=eqmap(i,j)/zod	!ratio, set a limit, & take the log
c			if(temp.le.0.001)temp=0.001
c		        temp=10.*alog10(temp)
c
			heliomap(ii,jj,n)=temp	!note that old values from previous cameras persist in here, if this value is zero
			nheliomap(ii,jj)=nheliomap(ii,jj)+1	!this counts contributions from all cameras, but note re-done below
c			
			jj=nint(10.*(rEloP+0.05))
			ii=nint(10.*(PA+180.05))	!note here PA is directly used, and not theta1
			if(ii.gt.3600)ii=3600		!happens occasionally when PA = 180.000	
			if(ii.lt.1)ii=1	
			heliomap1(ii,jj)=heliomap1(ii,jj)+temp	!PA vs rEloP
			nheliomap1(ii,jj)=nheliomap1(ii,jj)+1
			ii=jj
			jj=jjsv			!jjsv is Rlat
			if(iii.eq.+1)then
				ii=1801-ii	!"East" = Dusk
			else
				ii=1800+ii	!"West" = Dawn
			endif
			if(ii.gt.3600)ii=3600
			if(ii.lt.1)ii=1	
			heliomap2(ii,jj)=heliomap2(ii,jj)+temp	!rEloP vs Rlat
			nheliomap2(ii,jj)=nheliomap2(ii,jj)+1
		    endif			
94		endif	!zero map entries skip to here
	    enddo
	enddo
c	
95	enddo	!end loop on n
	enddo	!end loop on cameras
c
	print *,'Averaging'
	do j=1,1800
c
c	form averages
c	
		do i=1,3600
			isum=0
			sum=0.
			do n=nmin,nmax
			if(heliomap(i,j,n).ne.0.)then
				sum=sum+heliomap(i,j,n)	!this includes values from previous cameras, if the present camera was zero
				isum=isum+1
			endif
			enddo
			heliomapout(i,j)=0.
c			nheliomap(i,j)=isum	!use this when the just-above composite is wanted, over-rides above calculation
			if(isum.ge.1)heliomapout(i,j)=sum/float(isum)
c
c	This section gets the median in place of straight average above
c
c			isum=0
c			heliomapout(i,j)=0.
c			do k=1,kmax	
c				array(k)=1.
c				if(heliomap(i,j,k).ne.0.)then
c					array(k)=heliomap(i,j,k)
c					isum=isum+1
c				endif
c			enddo
c			call indexx(kmax,array,indx)
c			if(isum.ge.1)heliomapout(i,j)=array(indx(kmax+1-isum))
c			if(isum.ge.3)then		!3 or more, use median
c				karg=kmax-isum/2
c				heliomapout(i,j)=array(indx(karg))
c			elseif(isum.eq.2)then		!2, use average (rare)
c				heliomapout(i,j)=(array(indx(kmax))+array(indx(kmax-1)))/2.
c			elseif(isum.eq.1)then		!just one, use it
c				heliomapout(i,j)=array(indx(kmax))
c			endif
c
			if(nheliomap1(i,j).ge.1)heliomap1(i,j)=heliomap1(i,j)/float(nheliomap1(i,j))
			if(nheliomap2(i,j).ge.1)heliomap2(i,j)=heliomap2(i,j)/float(nheliomap2(i,j))
		enddo
	enddo
c
c	fill gaps, put in neighbors average
c	
99	continue
	print *,'Filling gaps'
	do j=2,1799
		do i=2,3599
			if(nheliomap(i,j).eq.0)then
				isum=0
				sum=0.
				do ii=i-1,i+1
				do jj=j-1,j+1
					if(nheliomap(ii,jj).ne.0)then
						isum=isum+1
						sum=sum+heliomapout(ii,jj)
					endif
				enddo
				enddo
				if(isum.gt.0)heliomapout(i,j)=sum/float(isum)
			endif
			if(nheliomap1(i,j).eq.0)then
				isum=0
				sum=0.
				do ii=i-1,i+1
				do jj=j-1,j+1
					if(nheliomap1(ii,jj).ne.0)then
						isum=isum+1
						sum=sum+heliomap1(ii,jj)
					endif
				enddo
				enddo
				if(isum.gt.0)heliomap1(i,j)=sum/float(isum)
			endif
			if(nheliomap2(i,j).eq.0)then
				isum=0
				sum=0.
				do ii=i-1,i+1
				do jj=j-1,j+1
					if(nheliomap2(ii,jj).ne.0)then
						isum=isum+1
						sum=sum+heliomap2(ii,jj)
					endif
				enddo
				enddo
				if(isum.gt.0)heliomap2(i,j)=sum/float(isum)
			endif
		enddo
	enddo
c
c	do n=1,29	
c	sum=0.
c	isum=0
c	do j=851,950
c	do i=1751,1850
c		if(heliomap(i,j,n).ne.0.)then
c			sum=sum+heliomap(i,j,n)
c			isum=isum+1
c		endif
c	enddo
c	enddo
c	print 30,n,dis(n),sum,isum
c	write(30,30)n,dis(n),sum,isum,sum/float(isum)
c30	format(i4,f10.5,f15.0,i8,f10.5)	
c	enddo
c
c	write output files
c
	if(mmap.ne.0)print *,'Writing output files'
	if(mmap.eq.1.or.mmap.eq.4)then
		open(12,file=outfile)
		if(inumb.eq.1)open(15,file=outfilen)
		write(12,12)
		if(inumb.eq.1)write(15,15)
c12		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'-.4 .4')		!for log plots of ratio
12		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'-10 10')		!for difference plots, ADUs
c12		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'-100 100')		!for difference plots, ADUs
15		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'0 30')		!for difference plots, numbers
		do j=1,1800		
		write(12,110)(heliomapout(i,j),i=1,3600)
		if(inumb.eq.1)write(15,115)(nheliomap(i,j),i=1,3600)
115		format(80i4)		
		enddo
		close(12)
		if(inumb.eq.1)close(15)
	endif
	if(mmap.eq.2.or.mmap.eq.4)then
		open(13,file=outfile1)
		write(13,13)
c13		format('DSAA',/,'3600 1800',/,'-180 180',/,'0 180',/,'-.4 .4')
13		format('DSAA',/,'3600 1800',/,'-180 180',/,'0 180',/,'-10 10')
c13		format('DSAA',/,'3600 1800',/,'-180 180',/,'0 180',/,'-100 100')
		do j=1,1800		
		write(13,110)(heliomap1(i,j),i=1,3600)
		enddo
		close(13)
	endif
	if(mmap.eq.3.or.mmap.eq.4)then
		open(14,file=outfile2)
		write(14,14)
c14		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'-.4 .4')
14		format('DSAA',/,'3600 1800',/,'-180 180',/,'-90 90',/,'-10 10')
		do j=1,1800		
		write(14,110)(heliomap2(i,j),i=1,3600)
		enddo
		close(14)
	endif
	stop 'done'
	end	
c*********************************************************************************	
	subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT)
	implicit none
c INPUTS:
c	ALFA	real	phase angle of the pole Z(new) in old coordinate system
c	BETA	real	polar angle of the pole Z(new) in old coordinate system
c	GAMMA	real	phase angle of X(new) in coordinate system (2)
c	PHI	real	phase angle in old coordinate system
c	RLAT	real	latitude in old coordinate system (-90<=RLAT<=90)
c OUTPUTS:
c	PHI	real	phase angle in new coordinate system (0<=PHI<360)
c	RLAT	real	latitude in new coordinate system
            real        ALFA
            real        BETA
            real        GAMMA
            real        PHI
            real        RLAT

        double precision        SINB
        double precision        COSB
        double precision        SINL
        double precision        COSL
        double precision        SINP
        double precision        COSP
        double precision        COST

        ! Needed for g77 when using goniod.f
        ! What happens if these are available as intrinsic functions??

        double precision        dsind
        double precision        dcosd
        double precision        dacosd
        double precision        datan2d

        if (abs(RLAT) .gt. 90.0)print *,'Bummer: RLAT = ',RLAT 	!call Say('rotate','E','invalid','latitude')

        PHI  = PHI-ALFA                 ! Rotation (1)

        SINB = dsind(dble(BETA))        ! Rotation (2)
        COSB = dcosd(dble(BETA))
        SINL = dsind(dble(RLAT))
        COSL = dcosd(dble(RLAT))
        SINP = dsind(dble(PHI ))
        COSP = dcosd(dble(PHI ))

        COST = SINL*COSB+COSL*SINB*COSP ! cos(polar angle)

        if (dabs(COST) .ne. 1.0d0) then	! changed abs --> dabs from Paul's, probably doesn't matter...
            SINP = COSL*SINP
            COSP = COSL*COSB*COSP-SINL*SINB

        else                            !Point on new Z-axis; PHI not defined

            !-------
            ! The following statement is inserted for the benefit of the
            ! synoptic map drawing program.
            ! If BETA=0 then SINB=0, COSB=+/-1, i.e. the z-axis stays the same
            ! or is flipped by 180 degrees. The input phase angle for points on
            ! the z-axis will remain the same or is changed to 180-(phase angle)
            ! (just as for points not on the z-axis).

            if (SINB .eq. 0.0d0) COSP = COSB*COSP

        end if

        PHI  = datan2d(SINP,COSP)
        PHI  = PHI-GAMMA                ! Rotation (3)
        PHI  = amod(amod(PHI,360.0)+360.0,360.0)!0<=PHI<360
        RLAT = 90.0d0-dacosd(COST)        ! acosd returns between 0 and 180 deg

        return
        end
c*********************************************************************************	
        subroutine POINT_ON_LOS(XLNG,XLAT,RP,ELOLD,ELNEW,iEorW)
        implicit none
cINPUTS:
c	(either topocentric or heliocentric)
c
c	XLNG	real	longitude (0<=XLNG<=360) and ...
c	XLAT	real	.. latitude (-90<=XLAT<=90) of line of sight (degrees)
c	RP	real	radial distance to point P
c			(in units of the observer-Sun distance)
c OUTPUTS:
c	(either heliocentric or topocentric)
c
c	XLNG	real	longitude (0<=XLNG<360) and ...
c	XLAT	real	.. latitude (-90<=XLAT<=90) of point P (degrees)
c	RP	real	distance of point P in units of observer-Sun distance
c	ELOLD	real	elongation in the old coordinate system (deg); 0<=ELOLD<=180
c	ELNEW	real	elongation in the new coordinate system 0<=ELNEW<=180
c	iEorW	integer	+1 or -1; indicates on which side (East or West) of
c			the Sun-Earth line the point is located
            real        XLNG
            real        XLAT
            real        RP
            real        ELOLD
            real        ELNEW
            integer     iEorW

        double precision        DLNG
        double precision        DLAT
        double precision        DRP
        double precision        R
        double precision        cosE
        double precision        sinP
        double precision        cosP

        ! Needed for g77 when using goniod.f
        ! What happens if these are available as intrinsic functions??

        double precision        dsind
        double precision        dcosd
        double precision        dasind
        double precision        dacosd
        double precision        datan2d

        DLNG = dble(XLNG)               ! Convert to double
        DLAT = dble(XLAT)
        DRP  = dble(RP  )
        if (RP .lt. 0.0) then             ! Negative RP:
            DLNG = 180.0d0+DLNG           ! .. change to opposite direction
            DLAT = -DLAT
            DRP  = -DRP
        end if

        cosE    = DLNG
        if (dabs(cosE) .gt. 180.0d0) cosE = cosE-dsign(1.0d0,cosE)*360.0d0
        iEorW   = dnint( dsign(1.0d0,cosE) )

        cosE    = dcosd(DLNG)*dcosd(DLAT)! Cos(old elongation)
        ELOLD   = dacosd(cosE)          ! Old elongation in [0,180]
        
        R       = 1d0+DRP*(DRP-2d0*cosE)! Cosine rule: (dist Sun-P)**2
        if (R .gt. 0d0) then
            R           = dsqrt(R)                      ! Dist Sun-P
            RP          = R                             ! Real*4 output
            sinP        = dsind(DLAT)*DRP/R             ! Sin(Latitude)
            sinP        = dmin1(1d0,dmax1(-1d0,sinP))
            XLAT        = sngl(dasind(sinP))            ! Latitude [-90,90]
            cosP        = 1d0-DRP*cosE
            cosP        = dmin1(R,dmax1(-R,cosP))       ! prop. cos(XLNG)
            ELNEW       = dacosd(cosP/R)                ! New elongation
            sinP        = -DRP*dcosd(DLAT)*dsind(DLNG)  ! prop. sin(XLNG)
            XLNG        = sngl(datan2d(sinP,cosP))
            XLNG        = amod(amod(XLNG,360.0)+360.0,360.0)! 0<=XLNG<360

        else
            RP          = 0.0           ! (just to prevent 'divide by zero'
            XLAT        = 0.0           ! condition)
            XLNG        = 0.0
            ELOLD       = 0.0
            ELNEW       = 0.0

        end if

        return
        end
c*********************************************************************************	
        subroutine ECLIPTIC_EQUATOR(ID,iYr,Doy,PHI,RLAT)
        implicit none
cINPUTS:
c	ID	integer	ID=0: ecliptic     ---> equatorial
c			ID=1: equatorial   ---> ecliptic
c	iYr	integer	year of current date
c	Doy	real	day of year, including a fraction for the time of day
c	PHI	real	ecliptic longitude for equinox of current date
c	RLAT	real	ecliptic latitude for equinox of current date
c OUTPUTS:
c	PHI	real	right ascension
c	RLAT	real	declination
            integer     ID
            integer     iYr
            real        Doy
            real        PHI,BETA	!BETA added to this list to satisfy implicit none
            real        RLAT

        double precision        EPS
        double precision        JD
        double precision        JEpoch

        call Julian(0,iYr,Doy,JD,JEpoch)        !Get date in Julian years
        EPS = 23.439291d0-0.00013004d0*(JEpoch-2000.0d0) !Angle (ecliptic-equator)

        BETA = sngl(EPS)                        ! .. to equatorial coord
        if (ID .eq. 1) BETA = -BETA             ! .. to ecliptic coord
        call rotate(90.0,BETA,-90.0,PHI,RLAT)

        return
        end
c*********************************************************************************	
        subroutine Julian(IDIN,iYr,Doy,JDio,JEpoch)
        implicit none
c INPUTS:
c	IDIN	integer		ID=0  date    --> JD, Julian epoch, Besselian epoch
c				=1  JD      --> Date, Julian epoch, Besselian epoch
c				=2  Julian epoch    --> Date, JD, Besselian epoch
c		(not active)	=3  Besselian epoch --> Date, JD, Julian epoch
c				Add 10: use modified Julian days (=JD-2400000.5)
c				Add 20: use days relative to Jan. 1, 2000 (=JD-2451544.5)
c  INPUTS/OUTPUTS: (depending on value of ID)
c	iYr	integer			year; the year xxxBC should be entered as -xxx+1.
c	Doy	real			day of year (including fraction for the time of day).
c	JD	double precision	Julian day
c	JEpoch	double precision	Julian epoch = time in Julian years
c*******data type specifications are here added on to satisfy implicit none*******
            integer             IDIN,ID,I,J,Laps,Leap,iDoy,JD0,I365,I1582
            integer             iYr
            real                Doy
            double precision    JDio
            double precision    JEpoch

        parameter               (I365  =  365)
        parameter               (I1582 = -418)  !1582-2000
        double precision        JD,Frac
        double precision        BEpoch
        double precision        JD2000

        JD2000 = 2451545        ! Jan 1.5, 2000 (noon)
        if (IDIN .ge. 10) JD2000 = 51544.5d0
        if (IDIN .ge. 20) JD2000 = 0.5d0

        ID = mod(IDIN,10)       ! ID=0,1,2,3

        if (ID .eq. 0) then     ! Date ---> Julian day, etc.
            Laps = iYr-2000     ! Whole years from 1 Jan 1.0d 2000 to idem of iYr
            I = Laps
            if (Laps .gt. 0) I = Laps-1
            J = max(I1582,I)
            Leap = I/4-J/100+J/400
            if (Laps .gt. 0) Leap = Leap+1      ! # leap years
            JD = I365*Laps+Leap+dble(Doy-1.5)   ! Rel. to 2000 Jan. 1.5 (noon)
            if (iYr .le. 1582) JD = JD+10

            JEpoch = 2000+JD/365.25d0
            BEpoch = 1900+(JD+36524.68648d0)/365.242198781d0

        else                            ! Julian day, etc. ---> Date
            if (ID .eq. 2) then         ! Julian epoch
                JD = JEpoch-2000
                BEpoch = 2000.0012775d0+1.000021359d0*JD
                JD = JD*365.25d0
            else if (ID .eq. 3) then
                JD = BEpoch-1900
                JEpoch = 1900.000858d0+0.999978642d0*JD
                JD = -36524.68648d0+JD*365.242198781d0
            else                        ! Assume ID=1
                JD = JDio-JD2000
                JEpoch = 2000+JD/365.25d0
                BEpoch = 1900+(JD+36524.68648d0)/365.242198781d0
            end if

            JD = JD+0.5d0               ! Shift to previous midnight
            JD0 = int(JD)               ! Whole days between previous midnight ...
            if (JD .lt. JD0) JD0 = JD0-1! .. and midnight of 1 Jan 2000
            Frac = JD-JD0               ! Day fraction from previous midnight
            Laps = 0
            iDoy = JD0
            do while (iabs(iDoy) .gt. 365)
                Laps = Laps+iDoy/366    ! Underestimate number of full years in iDoy
                I = Laps
                if (Laps .gt. 0) I = Laps-1
                J = max(I1582,I)
                Leap = I/4-J/100+J/400
                if (Laps .gt. 0) Leap = Leap+1  ! # leap years
                iDoy = JD0-Laps*I365-Leap
                if (Laps .le. I1582) iDoy = iDoy-10
            end do

            !-----
            ! Leave do while loop with -365<=iDoy<=365. Doy <= 0 happens for dates earlier
            ! than 2000 Jan. 1.0; Doy >=0 for dates later than 2000 Jan 1.0.

            iYr = 2000+Laps
            if (iDoy .lt. 0) then
                iYr = iYr-1
                if (iYr .eq. 1582) then
                    iDoy = iDoy+355     ! 355: # days in 1582
                    if (iDoy .lt. 0) iYr = iYr-1
                end if
            end if

            Leap = 0                    ! Suppose it's not a leap year
            if (iYr/4*4 .eq. iYr .and. (iYr .lt. 1582 .or. !It's a leap year
     &          iYr/100*100 .ne. iYr .or. iYr/400*400 .eq. iYr)) Leap = 1
            if (iDoy .lt. 0) iDoy = 365+Leap+iDoy       ! Now iDoy>=0

            if (iDoy .eq. 365+Leap) then
                iDoy = 0
                iYr = iYr+1
            end if

            Doy = iDoy+1+Frac
            JD = JD-0.5D0

        end if

        JDio = JD2000+JD

        return
        end
c*********************************************************************************	
        subroutine SunNewcomb(ID,iYr,Doy,rLng,rLat,rDis)
        implicit none
cINPUTS:
c	ID	integer		0 = return two-body solution
c				1 = take planetary perturbations into account
c	iYr	integer		year
c	Doy	real		doy of year, including fraction for time of day
c OUTPUTS:
c	rLng	double precision	ecliptic longitude in degrees [0,360)
c	rLat	double precision	ecliptic latitude  in arcsec
c	rDis	double precision	Sun-Earth distance in AU
            integer             ID
            integer             iYr
            real                Doy
            double precision    rLng
            double precision    rLat
            double precision    rDis

        double precision        JD
        double precision        T
        double precision        DLP
        double precision        G
        double precision        G2
        double precision        G4
        double precision        G5
        double precision        G6
        double precision        D
        double precision        A
        double precision        U
        double precision        DLO
        double precision        SLO
        double precision        DL
        double precision        R1

        ! Needed for g77 when using goniod.f
        ! What happens if these are available as intrinsic functions??

        double precision        dsind
        double precision        dcosd

        call Julian(0,iYr,Doy,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        
c*********************************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)	!from Numerical Receipes...
      implicit none
      integer*4 N
      integer*4 INDX(N),I,J,L,IR,INDXT
      real*4 ARRIN(N),Q
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END