C+
C NAME:
C	Extractdn
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.
C CALLING SEQUENCE:
C	call Extractdn(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,FALLOFF,
C    &		        VV,DD,XCshift,DVfact,DDfact,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		beginning time of observations
C       XCtend          real*8		ending time of observations
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# heights
C	nT		integer		# times
C	FALLOFF		real		Density falloff with solar distance
C	VV(nLng,nLat,nT)real		velocity map at height RR
C	DD(nLng,nLat,nT)real		density map at height RR
C	XCshift(nLng,nLat,nMap,nT,3) real shift amount from array
C	DVfact(nLng,nLat,nMap,nT) real  velocity difference at height		
C	DDfact(nLng,nLat,nMap,nT) real  density difference at height		
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 Extractdn(bExtract,cPrefix,NTP,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,NinterD,FALLOFF,
     &		            VV,DD,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)
	logical		bExtract(NTP)
	character 	cPrefix(NTP)*2
	real		VV(nLng,nLat,nT),
     &			DD(nLng,nLat,nT),
     &			XCshift(nLng,nLat,nMap,nT,3),
     &			DDfact (nLng,nLat,nMap,nT),
     &			DVfact (nLng,nLat,nMap,nT),
     &			dXC(3)

	real*8		JDCar(nCar),JD,JEpoch,JDCntr
	
	real*8		XCtbeg,XCtend,XCtfull,X,Doy8,ObsTim,ObsTimN
	
	real*8		dT				! Time increment for extracted time series
	real*8		DoypY		/0.0d0/
	
	real 		ALng		/3.0/
	
	real		RVect(3),XCbe(2,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
	
	external	EARTH8

	include		'MAPCOORDINATES.H'
	
C	print *, ' Into extractd.f XCstrt, XCtbeg, XCtend =',XCstrt, XCtbeg, XCtend

	Bad = BadR4()
	
	dT = (XCtend - XCtbeg)/(dble((NinterD+1)*(nT-1)))	! extract interval length in days
	
	NTT = nint(20.0/dT)					! 2*20 =# times extracted (potentially covers +-20 days from the center of the rotation)

C	print *, 'Extractdn Interval length in days', dT
	
	XCcntr = XCstrt + 0.5			! Carrington variable at center of rotation
	I = XCcntr				! Carrington interval before Earth was at center of rotation
	JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*dble(XCcntr-I)
	do ITP=1,NTP
	    if (bExtract(ITP)) then
						! 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
		write (cFile,'(A,A,F8.3)') cPrefix(ITP),'_',NCoff+XCstrt
		if (.not.bGetLun(iU)) call Say('Extract','E','NOUNIT',
     &			'iGetLun unable to assign logical unit number')		
		open (iU,file=cFile,status='NEW')
		call OSGetDateAndTime(cDate)
		write (iU,'(4A)') '; ',cFile(:itrim(cFile)),' created on ',cDate

		XCinc = 0.0
		iYrL  = 0
		DoypY = 0.0d0 
		do I=-NTT,NTT
	
		  call Julian8(1,iYr,Doy8,dble(JD+I*dT),JEpoch)

		  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

		  if((Doy8+DoypY).ge.XCtbeg.and.(Doy8+DoypY).le.XCtend) then		  
		    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))
		    
C		    if(iYr.gt.iYrL.and.I.gt.-NTT) then					! Year change
C	    		Leap = 0							! Suppose it's not a leap year
C	    		if ((iYrL)/4*4 .eq. (iYrL) .and. ((iYrL) .lt. 1582 .or. 	! It's a leap year
C     &	        	iYrL/100*100 .ne. iYrL .or. iYrL/400*400 .eq. iYrL)) Leap = 1
C     	    		DoypY = dble(365+Leap)
C			iYrL = iYr
C     		    end if
C		    iYrL = iYr

		    ObsTimN = ObsTim + DoypY
		    ObsTimN = min(XCtend,max(XCtbeg,ObstimN))

C	print *, 'ObsTimN',ObsTimN

		    call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,XCshift,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXC)
     
C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN dXC ', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN, dXC(1), dXC(2), dXC(3)

C		    if (dXC(1) .eq. Bad) stop 'Extract: oops 1'
C		    if (dXC(2) .eq. Bad) stop 'Extract: oops 2'
C		    if (dXC(3) .eq. Bad) stop 'Extract: oops 3'
		    if (dXC(1).ne.Bad.and.dXC(2).ne.Bad.and.dXC(3).ne.Bad) then

		      call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,DDfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,DF)
		      call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nt,DVfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,VF)

		      ObsXCN = ObsXC   + dXC(1)
		      ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN))

		      ObsLa  = ObsLat  + dXC(2)

		      ObsTimN = ObsTim + dXC(3)
		      ObsTimN = min(XCtend,max(XCtbeg,ObstimN))

		      call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DD,1,
     &			ObsXCN,ObsLa,ObsTimN,ObsD)
C		      print *, 'DF ObsD', DF, ObsD
		      if (DF .ne. Bad .and. ObsD .ne. Bad) ObsD = ObsD*DF*((RR/ObsDis)**FALLOFF)
		      call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1,
     &			ObsXCN,ObsLa,ObsTimN,ObsV)

C		      print *, 'After = ', JD+I*DT, iYr, DF, VF, Doy8, ObsXCN, ObsLa, ObsTimN, ObsDis, ObsD, ObsV

C		      print *, 'VF ObsV', VF, ObsV


		      if (VF .ne. Bad .and. ObsV .ne. Bad) ObsV = ObsV*VF
		      ObsD = max(ObsD,-9999.999)			! Should fit output format
		      ObsV = max(ObsV,-9999.999)
		    else
		      print *, cPrefix(ITP),' is out of range of the 3D matrix' 
		      ObsD = -9999.999
		      ObsV = -9999.999
		    end if

		    write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
		  end if
		end do
		iU = iFreeLun(iU)
	    end if
	end do
	return
	end
	
	subroutine ExtractPositionn8(ITP,iYr,Doy8,RVect)
	
	real		RVect(3),VVect(5)
	
	real*8		dLngSun,dLatSun,dDisSun
	
	real*8 		Doy8
	
	if (ITP .eq. 3) then		! Earth
	    call SunNewcomb8(0,iYr,Doy8,dLngSun,dLatSun,dDisSun)
	    RVect(1) =  dLngSun+180
	    RVect(2) = -dLatSun/3600	! arcsec -> degrees
	    RVect(3) =  dDisSun
	else if (ITP .eq. 10) then	! Ulysses
	    Doy = sngl(Doy8)
	    call UlyssesOrbit(iYr,Doy,RVect,VVect)
	else if (ITP .eq. 11 .or.	! Helios 1
     &		 ITP .eq. 12) then	! Helios 2
     	    Doy = sngl(Doy8)
	    call HOSOrbit(ITP-10,iYr,Doy,RVect(3),RVect(1),VR,VT)
	    RVect(2) = 0
	else if (ITP .eq. 13) then	! STEREO A
	    Doy = sngl(Doy8)
	    call stereoaorbit(iYr,Doy,RVect,VVect)
	else if (ITP .eq. 14) then	! STEREO B
	    Doy = sngl(Doy8)
	    call stereoborbit(iYr,Doy,RVect,VVect)
	else
	    Doy = sngl(Doy8)
	    call PlanetOrbit(ITP,iYr,Doy,RVect,VVect)
	end if
					! Convert to heliographic coordinates
	call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,RVect(1),RVect(2))
	return
	end