C+
C NAME:
C	Extractd3d
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.  This is the same program as Extracttd.for
C 	except that it extracts an interpolated time series from the density
C	matrix directly (using a 4-D interpolative scheme rather than from
C	the 3D base density using XCshift and DVfact and DDfact.  The densities
C	and velocities are extracted at UT times.  
C CALLING SEQUENCE:
C	call Extractd3D(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,
C    &			nMap,nT,nTn,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar)
C INPUTS:
C	bExtract(12)	logical		.TRUE.: extract; .FALSE.: do not extract
C		1<=ITP<=9 Planets (1=Mercury, 2=Venus, 3=Earth, etc.)
C		10=Ulysses, 11=Helios 1; 12=Helios 2
C	RR		real		height of deconvolution
C	dRR		real		distance between heights
C	XCbe(2,nT)	real		start and ending Carrington variable of arrays
C       XCtbeg          real*8		start time of observations
C       XCtend          real*8		end time of observations
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# heights
C	nT		integer		# of times
C	nTn		integer		# of times of extracted matrix
C	V3DT(nLng,nLat,nMap,nTn*nT) real  3 D velocity over time
C	D3DT(nLng,nLat,nMap,nTn*nT) real  3 D density over time
C	
C	NCoff		integer
C	XCstrt		real		Start rotation value
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C OUTPUTS:
C FUNCTIONS/SUBROUTINES:
C 	SunNewcomb, ECLIP_EQUATOR,WR2DARR
C-
	subroutine Extractd3D(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,
     &		nMap,nT,nTn,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar)
	parameter	(NTP = 12)
	logical		bExtract(NTP)
	character 	cPrefix(NTP)*2	/'ME','VE','E3','MA','JU','SA','UR','NE','PL','UL','H1','H2'/
	real		V3DT(nLng,nLat,nMap,nTn*nT),
     &			D3DT(nLng,nLat,nMap,nTn*nT)
	real*8		JDCar(nCar),JD,JEpoch,JDCntr,Doy8
	real		RVect(3),XCbe(2,nTn*nT)
	equivalence	(ObsLng,RVect(1)),(ObsLat,RVect(2)),(ObsDis,RVect(3))
	character 	cFile*14,
     &			cStr*90	/';Year DOY   Hour Height    Lat     Long   Density  Velocity'/,
     &			cDate*23
	logical		bGetLun
	parameter	(NTT = 60)		! 2*NT+1=# times extracted (covers 30 days)
	real*8		dT		/0.25d0/! Time increment for extracted time series
	real*8		XCtbeg,XCtend,XCtfull,X,ObsTim,ObsTimN
	real*8		DoypY		/0.0d0/
	
	include		'MAPCOORDINATES.H'

	Bad = BadR4()
	ObsTimNL = 0.0
	
	print *, ' INTO EXTRACTD3D XCstrt, XCtbeg, XCtend =',XCstrt, XCtbeg, XCtend
	
	XCcntr = XCstrt+0.5			! Carrington variable at center of rotation
	I = XCcntr				! Time when Earth was at center of rotation
	JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*(XCcntr-I)
	
	do ITP=1,NTP
	    if (bExtract(ITP)) then		! Check if this object is on the list to extract a position
						! Heliographic pos of extraction point at JDcntr
		call Julian8(1,iYr,Doy8,JDCntr,JEpoch)	
		call ExtractPositionn8(ITP,iYr,Doy8,RVect)
		
		JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*Sun__SynodicP
		I  = JD
		JD = I+nint((JD-I)/dT)*dT						! JD at the "center" of the extract location
		if (.not. bGetLun(iU)) call Say('Extract','E','NOUNIT',
     &			'iGetLun unable to assign logical unit number')
		write (cFile,'(A,A,F8.3)') cPrefix(ITP),'_',NCoff+XCstrt
		open (iU,file=cFile,status='NEW')
		call OSGetDateAndTime(cDate)
		write (iU,'(4A)') '; ',cFile(:itrim(cFile)),' created on ',cDate
		write (iU,'(A)' ) cStr(:itrim(cStr))
		XCinc = 0.0
		iYrL  = 0
		do I=-NTT,NTT
		    call Julian8(1,iYr,Doy8,JD+I*dT,JEpoch)
		    call ExtractPositionn8(ITP,iYr,Doy8,RVect)
		    ObsD = Bad
		    ObsV = Bad
		    ObsXC = XMAP_OBS_POS(XCstrt,ObsLng)
		    if(I.gt.-NTT) then				! Other times
			dObs = ObsXC-ObsXCL
		    else
			dObs = 0.
		    end if
		    ObsXCL = ObsXC
		    if(abs(dObs).gt.0.5) XCinc = XCinc + 1.0
		    ObsXC  = ObsXC + XCinc
C	print *, 'ObsXC', ObsXC
		    ObsTim = Doy8
C	print *, 'ObsTim',ObsTim
C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXC', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXC

		    ObsXCN = ObsXC
		    ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN))
		    	    
		    if(iYr.gt.iYrL.and.I.gt.-NTT) then					! Year change
	    		Leap = 0							! Suppose it's not a leap year
	    		if ((iYrL)/4*4 .eq. (iYrL) .and. ((iYrL) .lt. 1582 .or. 	! It's a leap year
     &	        	iYrL/100*100 .ne. iYrL .or. iYrL/400*400 .eq. iYrL)) Leap = 1
     	    		DoypY = dble(365+Leap)
			iYrL = iYr
     		    end if	
		    iYrL = iYr

		    ObsTimN = ObsTim + DoypY
		    ObsTimN = min(XCtend,max(XCtbeg,ObsTimN))     ! This is the interval time, not the time to interpolate

		    call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			  nMap,nTn*nT,D3DT,1,ObsTimN,ObsXCN,ObsLat,ObsDis,DF)
		    call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			  nMap,nTn*nT,V3DT,1,ObsTimN,ObsXCN,ObsLat,ObsDis,VF)
		    ObsD = Bad
		    if (DF .ne. Bad) ObsD = DF/(ObsDis*ObsDis)
		    if(ObsTimN .eq. ObsTimNL) ObsD = Bad
		    ObsV = Bad
		    if (VF .ne. Bad) ObsV = VF
		    if(ObsTimN .eq. ObsTimNL) ObsV = Bad
		    ObsTimNL = ObsTimN
		    ObsD = max(ObsD,-9999.999)			! Should fit output format
		    ObsV = max(ObsV,-9999.999)
C		    print *, iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsXCN,ObsLat,ObsTimN,ObsLng,dXc,DF,ObsD,VF,ObsV
		    write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV

		end do
		iU = iFreeLun(iU)
	    end if
	end do
	return
	end