C+
C NAME:
C	MkPostd
C PURPOSE:
C	To determine line of sight positions and their projections at a given 
C	height above the Sun.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkPostd(XCbe,XCtbeg,XCtend,Idaybeg,RRS,dRR,DLOS,nLng,nLat,nMap,nT,
C    &		XCshift,NL,NS1,DIST,XLS,XCE,XP,XD,IYRS,DOYS8,
C    &		XCproj,XClos,XLproj,XLlos,RRlos,XCtpr,XCtim,LON,LAT,ITIM)
C INPUTS:
C	XCbe(2,nT)	real		Time interval map beginnings and endings
C	XCtbeg		real		Time interval beginning
C	XCtend		real		Time interval ending
C	Idaybeg		integer		Total integer number of days in beginning year associated with XCtbeg
C	RRS		real		Distance above Sun for deconvolution maps
C	dRR		real		Heights between individual maps (in AU)
C	DLOS		real		resolution along LOS (in units of DIST)
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# TS MAPS (height begins at 0 AU)
C	nT		integer		# time intervals
C	XCshift(nLng,nLat,nMap,nT,3)	real 	Shifts to map to height RR (set by MkXCShifts)
C	NL		integer		# G-level data points
C	NS1		integer		# segments of length DLOS along LOS + 1
C					(the extra element is used to store Point-P info)
C	DIST(NL)	real		Observer-Sun distance at source time
C	XLS(NL)		real		Ecliptic longitude of Sun
C	XCE(NL)		real		Carrington variable of Earth
C	XP(NL)		real		Longitude diff. relative to Earth-Sun
C	XD(NL)		real		Declination of Source
C	IYRS(NL)	integer		Year of source observation
C	DOYS8(NL)	real*8		Day of year (and fraction) of source observation
C OUTPUTS:
C	XCproj(NS1,NL)	real		Modified Carrington variables for points on LOS,
C					after projection to reference surface RR
C	XClos(NS1,NL)	real		Modified Carrington variables for points on LOS
C	XLproj(NS1,NL)	real		Heliographic latitude for L.O.S. points on the reference surface
C	XLlos(NS1,NL)	real		Heliographic latitude for points on LOS
C	RRlos(NS1,NL)	real		Height above Sun for points on L.O.S.
C	XCtpr(NS1,NL)	real*8		Time of observation on reference surface
C	XCtim(NS1,NL)	real*8		Time of observation on L.O.S.
C	LON  (NS1,NL)	integer		Int value of the heliographic longitude
C	LAT  (NS1,NL)	integer		Int value of the heliographic latitude
C	ITIM (NS1,NL)	integer		Int value of time on reference surface
C FUNCTIONS/SUBROUTINES:
C	POINT_ON_LOS, ECLIPTIC_HELIOGRAPHIC, XMAP_OBS_POS
C PROCEDURE:
C >	XClos and XCproj are on the same scale as XCE.
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkPostd(XCbe,XCtbeg,XCtend,Idaybeg,RRS,dRR,DLOS,nLng,nLat,nMap,nT,
     &		XCshift,NL,NS1,DIST,XLS,XCE,XP,XD,IYRS,DOYS8,
     &		XCproj,XClos,XLproj,XLlos,RRlos,XCtpr,XCtim,LON,LAT,ITIM)

	real		XCshift(nLng,nLat,nMap,nT,3),! Map of TS shifts to RR map height
     &			XCbe(2,nt),		! Beginning and ending times
     &			DIST(NL),		! Distance of Sun to Earth
     &			XLS(NL),		! Ecliptic longitude of Sun
     &			XCE(NL),		! UT Time at Earth (in DOY)
     &			XP(NL),			! Longitude diff. relative to Earth-Sun
     &			XD(NL),			! Declination of Source
     &			XCproj(NS1,NL),		! Projected Carr. rot. value of point on L.O.S. 
     &			XClos (NS1,NL),		! Point on L.O.S. Carr. rot. value
     &			XLproj(NS1,NL),		! Latitude on the reference surface  
     &			XLlos (NS1,NL),		! Heliographic lat. (in deg.) point on L.O.S. 
     &			RRlos (NS1,NL),		! Height above Sun of point on L.O.S.
     &			XCbeg,XCend,		! Beginning and ending Carrington variable values
     &			XC(3)			! Internal array
     
        real*8		DOYS8(NL)		! Day of year (1.0 at beginning) (and fraction) of L.O.S. obs. (for SMEI, double precision)
     
     	real*8		Doy8,			! Day of Year (1.0 at beginning) to use for ECLIPTIC_HELIOGRAPHIC8
     &			XCtstrt,
     &			Adaybeg
     
     	real*8		XCtbeg,XCtend,XCtfull,X,tXC,XCtinc,T1
	
	real*8		XCtpr (NS1,NL),		! Projected time on reference surface
     &			XCtim (NS1,NL)		! Time in 3D space along LOS

	integer		IYRS(NL),		! Year of source observation
     &			LON(NS1,NL),		! Nearest longitude index
     &			LAT(NS1,NL),		! Nearest latitude index
     &			ITIM(NS1,NL)		! Nearest time index
	parameter	(DLOS25 = 0.25)
	include		'MAPCOORDINATES.H'
	XCinc	= (XCbe(2,1)-XCbe(1,1))/(2.*nLng1) ! Increment beyond beginning and end longitude
	XCtinc	= (XCtend-XCtbeg)/(2.D0*dble(nT1)) ! Increment beyond end of time
	XCbe1 	= XCbe(1,1)			! Beginning longitude limit
	XCtstrt	= XCtbeg			! Beginning time in DOY
	XCfin	= XCbe(2,nT)			! Ending longitude limit

C	print *, 'XCbe1, XCfin, XCtstrt', XCbe1,XCfin,XCtstrt
	do I=1,NL
	    iYr    = IYRS(I)
	    Doy8   = DOYS8(I)			! Double precision DOY value of line-of-sight observation
	    ObsR   = DIST(I)			! Heliocentric distance observer
	    ObsLng = XLS(I)-180.		! Heliocentric ecliptic longitude Earth
	    ObsXC  = XCE(I)			! Modified Carrington variable Earth
C	    If(I.lt.10) print *, 'Inside mkpostd',IYRS(I),DOYS8(I),DIST(I),ObsLng,OBsxc
C	    tXC	   = DOYS(I)-XCtstrt		! Time of source observation (days since start) in single precision
	    tXC	   = DOYS8(I)			! Time of source observation (days since year beginning) in double precision
C	    if(XCtend.gt.367.00) then		! Handle case where one year to another is crossed
C	    	iYr = iYr - 1	    
C	    end if

	    Adaybeg = dble(Idaybeg)		! This is the number of days in the preceeding year
	    if(XCtend.gt.Adaybeg.and.DOYS8(I).lt.dble(180.)) then ! Handle case where one year to another is crossed
	        tXC  = DOYS8(I) + Adaybeg - 1.0d0
	        Doy8 = Doy8     + Adaybeg - 1.0d0
	        iYr = iYr - 1			! The source year time should be one greater, so subtract one year
	    end if
C	    If(I.lt.3) then
C	      print *, 'Inside mkpostd',IYRS(I),iYr,tXC,DOYS8(I),DIST(I),ObsLng,OBsXC,Adaybeg
C	      print *,XCinc,XCbe(2,1),XCbe(1,1),XCtinc,XCtend,Xctbeg,XCfin,Idaybeg
C	    end if
	    
	    do J=1,NS1				! Loop over all segments along LOS
						! Define point P along LOS:
		P1 = XP(I)			! Topocentric ecliptic longitude diff. with Earth-Sun line
		P2 = XD(I)			! Topocentric ecliptic latitude
		P3 = (J-.5)*DLOS		! Distance along LOS (units of ObsR)
						! Point-P distance along LOS (units of ObsR)
		if (J .eq. NS1) P3 = max(DLOS25,cosd(P2)*cosd(P1))
						! Convert to heliocentric coordinates
		T1 = tXC-dble(P3*ObsR*TCD)	! Time of L.O.S. segment obs.
		T1 = min(XCtend+XCtinc,max(XCtbeg-XCtinc,T1)) ! No T1 outside limits
		XCtim(J,I) = T1
		call POINT_ON_LOS(P1,P2,P3,ELO,ELOS,iEorW)
		P3 = P3*ObsR			! Convert to AU
		P1 = ObsLng+P1			! Heliocentric longitude of P
		call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,P1,P2)
						! Convert to heliographic coordinates
		P1 = XMAP_OBS_POS(ObsXC,P1)	! Modified Carrington variable of P
		P1 = min(XCfin+XCinc,max(XCbe1-XCinc,P1)) ! No P1 outside longitude limits
C		print *, ' '
C		print *, 'before Get4Dval - I J T1 P1 P2 P3', I,J, T1,P1,P2,P3
		call Get4Dval(3,XCbe,XCtbeg,XCtend,RRS,dRR,nLng,nLat,nMap,nT,
     &			XCshift,1,T1,P1,P2,P3,XC)
C		if(J.eq.20.and.I.eq.1000)print *, 'after Get4Dval - XC =',  XC(1), XC(2), XC(3)
		XClos(J,I) = P1			! Modified Carrington variable of P
		XLlos(J,I) = P2			! Heliographic latitude of P
		RRlos(J,I) = P3			! Heliocentric distance of P (AU)

		XCproj(J,I)= P1+XC(1)		! Projected modified Carrington variable
		XCproj(J,I)=min(XCfin+XCinc,max(XCbe1-XCinc,XCproj(J,I))) ! No XCproj outside longitude limits

		XLproj(J,I)= P2+XC(2)		! Projected heliographic latitude of P
		XLproj(J,I)=max(-90.,min(90.,XLproj(J,I))) ! No XLproj outside latitude limits	
		
		XCtpr(J,I) = T1+dble(XC(3))  	! Proj. time from beginning CL value		
		XCtpr(J,I) = min(XCtend+XCtinc,max(XCtbeg-XCtinc,XCtpr(J,I))) ! No XCtpr outside limits

		ITIM(J,I)= nint(XCtvar(XCtpr(J,I))) ! Nearest interval time index
C	if(ITIM(J,I).eq.3.and.J.eq.1) print *, ' '
C	if(ITIM(J,I).eq.3.and.J.eq.1) print *, ' ITIM', ITIM(J,I), J
C	if(ITIM(J,I).eq.3.and.J.eq.NS1) print *, ' ITIM', ITIM(J,I), J

C	print *, ITIM(J,I),XCtpr(J,I),T1,dble(XC(3)),T1+dble(XC(3))

		ITIM(J,I)= max(1,min(nT,ITIM(J,I))) ! No IT outside limits

		LAT(J,I) = nint(XLindx(XLproj (J,I))) ! Nearest latitude index
		LAT(J,I) = max(1,min(nLat,LAT(j,I)))  ! No LAT outside limits

		XCbeg	 = XCbe(1,ITIM(J,I)) 
		XCend	 = XCbe(2,ITIM(J,I))
		LON(J,I) = nint(XCvar(XCproj(J,I)))  ! Nearest longitude index
		LON(J,I) = max(1,min(nLng,LON(J,I))) ! No LON outside limits

C	if(ITIM(J,I).eq.3.and.J.eq.1) then
C	if(I.eq.10.and.J.eq.1) then
C	print *, 'At end of Mkpos.', J,I,XCproj(J,I),XClos(J,I),XLproj(J,I),XLlos(J,I),
C     &	RRlos(J,I),XCtpr(J,I),XCtim(J,I),LON(J,I),LAT(J,I),ITIM(J,I)
C	print *, 'At end of Mkpos.', XCbeg, XCend, XCtpr(J,I), 
C     &		XCtvar(XCtpr(J,I)), ITIM(J,I),XCproj(J,I),LON(J,I), LAT(J,I), J
C	end if
C	if(ITIM(J,I).eq.3.and.J.eq.NS1) then  
C	if(I.eq.10.and.J.eq.NS1) then
C	print *, 'At end of Mkpos.', XCbeg, XCend, XCtpr(J,I), 
C     &		XCtvar(XCtpr(J,I)), ITIM(J,I),XCproj(J,I),LON(J,I), LAT(J,I), J
C	end if
	    end do
	end do
	return
	end