How to use the jplephem.object_doppler function in jplephem

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):

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, \
				    N_Marsh, Marsh_alg, min_extract_col,\

gobs      = ephem.Observer()  
    gobs.name = h[0].header['TELESCOP']
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude)
    timeT = h[0].header['UTC-OBS'].split(':')
    if len(timeT[0]) == 1:
        gobs.date = h[0].header['DATE-OBS'][:10] + ' 0' + h[0].header['UTC-OBS']
    else:
        gobs.date = h[0].header['DATE-OBS'][:10] + ' ' + h[0].header['UTC-OBS']
    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:'

    if mode == 'so':
        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 ):

            sci_Ss = GLOBALutils.simple_extraction(data,c_all,ext_aperture,\

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 = 'Eso2.2'  
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude) 
    gobs.date = h[0].header['DATE-OBS'][:10] + ' ' + h[0].header['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
 
    ThAr_Ne_ref_m       = ThAr_Ne_ref
    ThAr_Ne_ref_dates_m = ThAr_Ne_ref_dates

    sorted_indices     = np.argsort( np.abs( np.array(ThAr_Ne_ref_dates_m) - mjd ) )
    sci_fits_ob        = dirout + fsim.split('/')[-1][:-4]+'spec.ob.fits.S'
    sci_fits_co        = dirout + fsim.split('/')[-1][:-4]+'spec.co.fits.S'
    sci_fits_ob_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple.ob.fits.S'
    sci_fits_co_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple.co.fits.S'

    if ( os.access(sci_fits_ob,os.F_OK) == False )        or ( os.access(sci_fits_co,os.F_OK) == False )  or \
    ( os.access(sci_fits_ob_simple,os.F_OK) == False ) or ( os.access(sci_fits_co_simple,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)

    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]) )
        sci_Ss = np.zeros( (nord,data.shape[1]) )

if len(Mmi)<2:
        Mmi = '0'+Mmi
    segs = (mins - int(Mmi))*60.
    if segs<10:
        Mse = '0'+str(segs)[:5]
    else:
        Mse = str(segs)[:6]

    gobs.date = str(DDATE[:4]) + '-' +  str(DDATE[4:6]) + '-' + str(DDATE[6:]) + ' ' +  Mho + ':' + Mmi +':' +Mse

    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 ) )
    sorted_indices_FP = np.argsort( np.abs( np.array(ThFP_ref_dates) - mjd ) )

    print '\t\tExtraction:'
    # optimally and simply extract spectra
    sci_fits_ob = dirout + fsim.split('/')[-1][:-8]+'spec.ob.fits.S'
    sci_fits_co = dirout + fsim.split('/')[-1][:-8]+'spec.co.fits.S'
    sci_fits_ob_simple = dirout + fsim.split('/')[-1][:-8]+'spec.simple.ob.fits.S'
    sci_fits_co_simple = dirout + fsim.split('/')[-1][:-8]+'spec.simple.co.fits.S'
    sci_fits_bac = dirout + fsim.split('/')[-1][:-8]+'spec.simple.bac.fits.S'
    if ( os.access(sci_fits_ob,os.F_OK) == False ) or \
        ( os.access(sci_fits_co,os.F_OK) == False ) or \

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 = 'VLT'
    gobs.lat  = rad(latitude)  # lat/long in decimal degrees  
    gobs.long = rad(longitude) 
    gobs.date = h[0].header['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

    sci_fits        = dirout + fsim.split('/')[-1][:-4]+'spec.pkl'
    sci_fits_simple = dirout + fsim.split('/')[-1][:-4]+'spec.simple.pkl'
    P_fits          = dirout + 'P_' + fsim.split('/')[-1][:-4]+'spec.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..."
        #"""
        print "\t\t\tweights for chip1..."
        P1 = GLOBALutils.obtain_P(data[:,:,0].T,c_all1,ext_aperture,ron1,gain1,NSigma_Marsh, S_Marsh, \
                                    N_Marsh, Marsh_alg, min_extract_cols1,max_extract_cols1, npools)
        print "\t\t\tweights for chip2..."

#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

	    data  = np.fliplr(data.T)

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']

        orders_offset = pdict['orders_offset']