How to use the jplephem.barycentric_object_track function in jplephem

# Moon Phase Calculations
    gobs = ephem.Observer()  
    gobs.name='Eso3.6'  
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude)
    DDATE     = h[0].header['DATE-OBS'][:10]
    HHOUR     = h[0].header['DATE-OBS'][11:]
    gobs.date = str(DDATE[:4]) + '-' +  str(DDATE[5:7]) + '-' + str(DDATE[8:]) + ' ' +  HHOUR[:2] + ':' + HHOUR[3:5] +':' + str(float(HHOUR[6:]) + halfcounts * TEXP )
    mephem    = ephem.Moon()
    mephem.compute(gobs)

    mephem = ephem.Moon()
    mephem.compute(gobs)
    Mcoo = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Mp = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Sp = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
    refvel = bcvel_baryc + moonvel
    print '\t\tRadial Velocity of sacttered moonlight:',refvel

    sorted_indices = np.argsort( np.abs( np.array(ThAr_ref_dates) - mjd ) )

    # optimally and simply extract spectra
    sci_fits_ob_B = dirout + fsim.split('/')[-1][:-4]+'spec.ob.B.fits.S'
    sci_fits_co_B = dirout + fsim.split('/')[-1][:-4]+'spec.co.B.fits.S'
    sci_fits_ob_R = dirout + fsim.split('/')[-1][:-4]+'spec.ob.R.fits.S'
    sci_fits_co_R = dirout + fsim.split('/')[-1][:-4]+'spec.co.R.fits.S'
    #force_sci_extract = True
    if ( os.access(sci_fits_ob_B,os.F_OK) == False ) or ( os.access(sci_fits_co_B,os.F_OK) == False ) \
    or ( os.access(sci_fits_ob_R,os.F_OK) == False ) or ( os.access(sci_fits_co_R,os.F_OK) == False ) \
    or (force_sci_extract):

print "\t\tBarycentric velocity:", bcvel_baryc

    res  = jplephem.pulse_delay(ra/15.0, dec, int(mjd), mjd%1, 1, 0.0)
    mbjd = mjd + res['delay'][0] / (3600.0 * 24.0)

    # Moon Phase Calculations
    gobs      = ephem.Observer()  
    gobs.name = 'Keck'
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude) 
    gobs.date = h[0].header['DATE-OBS'] + ' ' + h[0].header['UTC'].replace(':','-')
    mephem    = ephem.Moon()
    mephem.compute(gobs)
    Mcoo        = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Mp   = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Sp   = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
    refvel = bcvel_baryc + moonvel
    print '\t\tRadial Velocity of sacttered moonlight:',refvel
 
    sci_fits        = dirout + fsim.split('/')[-1][:-4]+'spec_'+str(int(chip))+'.fits'
    sci_fits_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple_'+str(int(chip))+'.fits'
    P_fits          = dirout + 'P_' + fsim.split('/')[-1][:-4]+'spec_'+str(int(chip))+'.fits'

    if ( os.access(sci_fits,os.F_OK) == False ) or ( os.access(sci_fits_simple,os.F_OK) == False ) or \
       ( force_sci_extract ):
	
        print "\t\tNo previous extraction or extraction forced for science file", fsim, "extracting..."

	P = GLOBALutils.obtain_P(data,c_new,ext_aperture,ronoise,\
                                    gain,NSigma_Marsh, S_Marsh, \

lbary_ltopo = 1.0 + res['frac'][0]
    bcvel_baryc = ( lbary_ltopo - 1.0 ) * 2.99792458E5
    print "\t\tBarycentric velocity:", bcvel_baryc
    res  = jplephem.pulse_delay(ra/15.0, dec, int(mjd), mjd%1, 1, 0.0)
    mbjd = mjd + res['delay'][0] / (3600.0 * 24.0)

    gobs      = ephem.Observer()  
    gobs.name = h[0].header['TELESCOP']
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude) 
    gobs.date = h[0].header['DATE'][:10] + ' ' + h[0].header['DATE'][11:]
    mephem    = ephem.Moon()
    mephem.compute(gobs)

    Mcoo = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Mp = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Sp = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
    refvel = bcvel_baryc + moonvel
    print '\t\tRadial Velocity of sacttered moonlight:',refvel
    moon_alts.update({fsim:mephem.alt})
    moon_ills.update({fsim:lunation})

    print '\t\tExtraction:'
    sci_fits        = dirout + fsim.split('/')[-1][:-4]+'spec.fits.S'
    sci_fits_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple.fits.S'
    if ( os.access(sci_fits,os.F_OK) == False ) or ( os.access(sci_fits_simple,os.F_OK) == False ) or \
       ( force_sci_extract ):

        print "\t\t\tNo previous extraction or extraction forced for science file", fsim, "extracting..."
        sci_S  = np.zeros( (nord,3,data.shape[1]) )

print "\t\tBarycentric velocity:", bcvel_baryc

    res  = jplephem.pulse_delay(ra/15.0, dec, int(mjd), mjd%1, 1, 0.0)
    mbjd = mjd + res['delay'][0] / (3600.0 * 24.0)

    # Moon Phase Calculations
    gobs      = ephem.Observer()  
    gobs.name = 'Keck'
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude) 
    gobs.date = h[0].header['DATE-OBS'] + ' ' + h[0].header['UTC'].replace(':','-')
    mephem    = ephem.Moon()
    mephem.compute(gobs)
    Mcoo        = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Mp   = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
    Sp   = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
    refvel = bcvel_baryc + moonvel
    print '\t\tRadial Velocity of sacttered moonlight:',refvel
 
    sci_fits        = dirout + fsim.split('/')[-1][:-4]+'spec_'+str(int(chip))+'.fits'
    sci_fits_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple_'+str(int(chip))+'.fits'
    P_fits          = dirout + 'P_' + fsim.split('/')[-1][:-4]+'spec_'+str(int(chip))+'.fits'

    if ( os.access(sci_fits,os.F_OK) == False ) or ( os.access(sci_fits_simple,os.F_OK) == False ) or \
       ( force_sci_extract ):
	
        print "\t\tNo previous extraction or extraction forced for science file", fsim, "extracting..."

	P = GLOBALutils.obtain_P(data,c_new,ext_aperture,ronoise,\

res         = jplephem.pulse_delay(RA/15.0, DEC, int(scmjd), scmjd%1, 1, 0.0)
        scmbjd      = scmjd + res['delay'][0] / (3600.0 * 24.0) #This is the modified barycentric julian day of the observation

        # set observatory info to retrive info about the moon
        gobs = ephem.Observer()
        gobs.name = 'VBT'
        gobs.lat  = rad(latitude)
        gobs.long = rad(longitude)
        #gobs.date = hd['UT-DATE'] + ' ' + hd['UT-TIME'].replace(':','_')
        gobs.date = hd['DATE-OBS'].replace('T',' ')

        mephem = ephem.Moon()
        mephem.compute(gobs)
        Mcoo = jplephem.object_track("Moon", int(scmjd), float(scmjd%1), 1, 0.0)
        Mp   = jplephem.barycentric_object_track("Moon", int(scmjd), float(scmjd%1), 1, 0.0)
        Sp   = jplephem.barycentric_object_track("Sun", int(scmjd), float(scmjd%1), 1, 0.0)
        res      = jplephem.object_doppler("Moon", int(scmjd), scmjd%1, 1, 0.0)
        lunation,moon_state,moonsep2,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,RA,DEC)
        refvel = bcvel_baryc + moonvel  #This is the velocity of the spectrum of the moon with the applied barycentric correction in the direction of the target. 

        print '\t\t\tBarycentric velocity:',refvel
        
        obj_fits        = dirout+nama+'.spec.fits.S'
        obj_fits_simple = dirout+nama+'.spec.simple.fits.S'
        obj_S           = pyfits.getdata(obj_fits)
        obj_Ss          = pyfits.getdata(obj_fits_simple)

        hth = pyfits.getheader(ThAr_ref[0])
        wavsol_pkl       = dirout+'ThAr_'+hth['DATE-OBS']+'.wavsolpars.pkl'
        pdict      = pickle.load(open(wavsol_pkl,'r'))
        global1     = pdict['p1']

gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
	    gobs.long = rad(longitude)

	    #date = hd['DATE-OBS']
	    #date = datetime.datetime(int(date[:4]),int(date[5:7]),int(date[8:10]),int(date[11:13]),int(date[14:16]),int(date[17:19]))
	    #new_date = date
	    #OJO aquiiiiii
	    #print 'Warning!!! adding 5 hrs to comute MJD due to problem in header! CHECK in future!!'
	    #new_date = date + datetime.timedelta(hours=5)

	    #gobs.date = new_date.strftime('%Y-%m-%d %H:%M:%S')
	    gobs.date = hd['DATE-OBS'][:10] + ' ' + hd['DATE-OBS'][11:]
	    mephem = ephem.Moon()
	    mephem.compute(gobs)
	    Mcoo = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	    Mp   = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	    Sp   = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
	    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
	    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
	    refvel = bcvel_baryc + moonvel
	    print '\t\tRadial Velocity of sacttered moonlight:',refvel

	    if os.access(dirout+'Dark_'+str(exptime)+'.fits',os.F_OK):
		DARK = pyfits.getdata(dirout+'Dark_'+str(exptime)+'.fits')
	    else:
		DARK = np.zeros(MasterBias[0].data.shape)

	    hdat  = pyfits.open( fsim )
	    data  = pyfits.getdata( fsim )
	    if data.shape == (2048, 2064):
		data = data.T
	    data = data - MasterBias[0].data #- DARK

gobs.long = rad(longitude)

	    #date = hd['DATE-OBS']
	    #date = datetime.datetime(int(date[:4]),int(date[5:7]),int(date[8:10]),int(date[11:13]),int(date[14:16]),int(date[17:19]))
	    #new_date = date
	    #OJO aquiiiiii
	    #print 'Warning!!! adding 5 hrs to comute MJD due to problem in header! CHECK in future!!'
	    #new_date = date + datetime.timedelta(hours=5)

	    #gobs.date = new_date.strftime('%Y-%m-%d %H:%M:%S')
	    gobs.date = hd['DATE-OBS'][:10] + ' ' + hd['DATE-OBS'][11:]
	    mephem = ephem.Moon()
	    mephem.compute(gobs)
	    Mcoo = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	    Mp   = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	    Sp   = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
	    res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
	    lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
	    refvel = bcvel_baryc + moonvel
	    print '\t\tRadial Velocity of sacttered moonlight:',refvel

	    if os.access(dirout+'Dark_'+str(exptime)+'.fits',os.F_OK):
		DARK = pyfits.getdata(dirout+'Dark_'+str(exptime)+'.fits')
	    else:
		DARK = np.zeros(MasterBias[0].data.shape)

	    hdat  = pyfits.open( fsim )
	    data  = pyfits.getdata( fsim )
	    if data.shape == (2048, 2064):
		data = data.T
	    data = data - MasterBias[0].data #- DARK

bcvel_baryc = 0.
		mbjd = mjd

	# Moon Phase Calculations
	gobs      = ephem.Observer()  
	gobs.name = 'ODUC'  
	gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
	gobs.long = rad(longitude)

	gobs.date = hd['DATE-OBS'].replace('T',' ')

	mephem = ephem.Moon()
	mephem.compute(gobs)
	Mcoo = jplephem.object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	Mp   = jplephem.barycentric_object_track("Moon", int(mjd), float(mjd%1), 1, 0.0)
	Sp   = jplephem.barycentric_object_track("Sun", int(mjd), float(mjd%1), 1, 0.0)
	res  = jplephem.object_doppler("Moon", int(mjd), mjd%1, 1, 0.0)
	lunation,moon_state,moonsep,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,ra,dec)
	refvel = bcvel_baryc + moonvel
	print '\t\tRadial Velocity of sacttered moonlight:',refvel


	#print utils.search_name(obj)
	nama = pucherosutils.search_name(obj)+'_'+hd['DATE-OBS'][:10]+'_'+hd['DATE-OBS'][11:13]+'-'+hd['DATE-OBS'][14:16]+'-'+hd['DATE-OBS'][17:]		
	obj_fits = dirout+nama+'.spec.fits.S'
	bkg_obj_fits = dirout+'Bkg_'+nama+'.fits'
	obj_fits_simple = dirout+nama+'.spec.simple.fits.S'
	

	if ( os.access(obj_fits,os.F_OK) == False )  or \
	   ( os.access(obj_fits_simple,os.F_OK) == False ) or \
	   (force_sci_extract):

bcvel_baryc = ( lbary_ltopo - 1.0 ) * 2.99792458E5  #This in the barycentric velocity
        res         = jplephem.pulse_delay(RA/15.0, DEC, int(scmjd), scmjd%1, 1, 0.0)
        scmbjd      = scmjd + res['delay'][0] / (3600.0 * 24.0) #This is the modified barycentric julian day of the observation

        # set observatory info to retrive info about the moon
        gobs = ephem.Observer()
        gobs.name = 'VBT'
        gobs.lat  = rad(latitude)
        gobs.long = rad(longitude)
        #gobs.date = hd['UT-DATE'] + ' ' + hd['UT-TIME'].replace(':','_')
        gobs.date = hd['DATE-OBS'].replace('T',' ')

        mephem = ephem.Moon()
        mephem.compute(gobs)
        Mcoo = jplephem.object_track("Moon", int(scmjd), float(scmjd%1), 1, 0.0)
        Mp   = jplephem.barycentric_object_track("Moon", int(scmjd), float(scmjd%1), 1, 0.0)
        Sp   = jplephem.barycentric_object_track("Sun", int(scmjd), float(scmjd%1), 1, 0.0)
        res      = jplephem.object_doppler("Moon", int(scmjd), scmjd%1, 1, 0.0)
        lunation,moon_state,moonsep2,moonvel = GLOBALutils.get_lunar_props(ephem,gobs,Mcoo,Mp,Sp,res,RA,DEC)
        refvel = bcvel_baryc + moonvel  #This is the velocity of the spectrum of the moon with the applied barycentric correction in the direction of the target. 

        print '\t\t\tBarycentric velocity:',refvel
        
        obj_fits        = dirout+nama+'.spec.fits.S'
        obj_fits_simple = dirout+nama+'.spec.simple.fits.S'
        obj_S           = pyfits.getdata(obj_fits)
        obj_Ss          = pyfits.getdata(obj_fits_simple)

        hth = pyfits.getheader(ThAr_ref[0])
        wavsol_pkl       = dirout+'ThAr_'+hth['DATE-OBS']+'.wavsolpars.pkl'
        pdict      = pickle.load(open(wavsol_pkl,'r'))
        global1     = pdict['p1']