mergefermi.py

# Python version of Micheal's fermimerge IDL program... hopefully. 

# Import as needed:

import subprocess
import time
import os
import pdb




# start timer:

start_time = time.time()

# functions for set up:
                                    
def check_for(x):
    f = os.getcwd().split('\\')[-1]
    x1 = '/' + x
    if not os.path.exists(f + x1):
        return True


def make_needed(x): 
    if check_for(x) == True:
        f = os.getcwd().split('\\')[-1]
        y = '.txt'
        if y in x:
            open(x, 'w').close()
            print ('made - ' + x)
        else:    
            os.makedirs(f +'/' + x)
            print ('made new directory - ' + x)

def merge(marker, new):
    if check_for(new) == True:
        g = 'ls *'
        i = '* > '
        h = g + marker + i  + new 
        subprocess.call(h, shell=True)
        print ('made ' + new)
    
def link(marker, new):
    if check_for(new) == True:
        g = 'ln -s *'
        i = '* > '
        h = g + marker + i  + new 
        subprocess.call(h, shell=True)
        print ('made ' + new)
        
def days_to_secs(days):
    secs = float(days*24*3600)
    return secs

def dph(orbitalperioddays): 
    dph = round(float(days_to_secs(orbitalperioddays)/nbins))
    return dph
    # calculates phase bin size
    
def ref_phase(days):
    ref1 = days_to_secs(float(days) - 51910.0)
    return ref1   
    #converts MJD to Fermi MET in secs

                        
# default values:


evclass = 2 # event class - do not change without good reason

zmax = 100 # max zenith angle - do not change without good reason

roi = 20.0 # region of interest - change as needed but have been using 20.0

emin = 100 # min energy - change as needed

emax = 300000 # max energy - change as needed

nbins = 10 # default number of bins per orbital phase

pshift = 0.0 #default phase to begin calculation

chatter = 1

# generic set up, input required later:

ra = 0.0  # degrees
dec = 0.0  # degrees
period = 1.0 # days
fileroot = '_sourcename' 
meti = 239557417 # value if using "START" for data acquisition
ref1 = meti #starting point for first loop of phase separation - will overwrite as 
#fermimerge iterates
metf = 255398400

# generic set up, runs automatically: 

pbs = dph(period)

period = float(days_to_secs(period))

make_needed('lightcurves')

merge('PH', 'ph.txt')

link('SC', 'sc.fits')

lc_dir = os.getcwd().split('\\')[-1] + '/lightcurves'





# functions that perform the work:


def dir_in(name, loc):
    os.chdir(loc) 
    make_needed(name)
    os.chdir('..')
    # creates new directory in a specified location (rather than in the cwd)



def gtselect(loc, i, tstart,tstop, wd):
    # pdb.set_trace()
    txt1 = 'gtselect evclass=' + str(evclass) + ' infile=@ph.txt  outfile=' 
    txt2 =  loc + 'events_ii_' + i + '.fits ra=INDEF dec=INDEF rad=INDEF tmin='
    txt3 =  str(tstart) + ' tmax=' + str(tstop) + ' emin=' + str(emin)  
    txt4 = ' emax=' + str(emax) + ' zmax=' + str(zmax) + ' chatter = 2'
    txt = txt1 + txt2 + txt3 + txt4
    print txt
    chd = 'cd' + wd
    #subprocess.call('ur_setup', shell = True)
    subprocess.call(chd, shell = True)
    subprocess.call(txt, shell = True)
    print ('made selection file')


def fermimerge(x):
   # p loop determines current phase and does set up for the next loop's work:
   p = x 
   if p < nbins + 1: 
      phase = p * pbs/period + pshift
      phase1 = (p + 1) * pbs/period + pshift
      print ('phase bin ' + str(phase) + '---' + str(phase1)) 
      phase = str(phase)
      bin_dir = 'ph_' + phase
      dir_in(bin_dir, lc_dir)
      tmploc = os.getcwd().split('\\')[-1] + '/lightcurves/' + bin_dir
      dir_in('tmp',tmploc)
      os.chdir('..')
      dfr = lc_dir + bin_dir + fileroot
      tstart = ref1
      
      # i loop defines beginning and end of phase bins and cycles 
      i = 00
      
      # first check that times are valid and adjust if needed:
      
      if tstart < metf + 1:
          tstart = ref1 + period*i + pbs*p + pshift*period 
          tstop = tstart + pbs
          select = 'no'
          if tstart >= meti and tstart <= metf:
              select = 'yes1'
          if tstop >= meti and tstop <= metf:
              select = 'yes2'
          if select == 'yes1' or 'yes2':
              if tstart < meti:
                  tstart = meti
              if tstop > metf:
                  tstop = metf                  
              print (tstart, tstop, tstop-tstart)
          
          # formatting for gtselect output file: 
          
          tmp = tmploc + '/tmp/' 
          if i < 10:
              i = '00' + str(i)
          if i >= 10 and i < 100:
              i = '0' + str(i)
          if i >= 100 and i < 1000:
              i = str(i)
          loc = os.getcwd()          
          gtselect(tmp, i, tstart, tstop, loc)    
                    
      

    
    
    
    

# and finally... : 

#ra = 0.0  # degrees
#dec = 0.0  # degrees
#period = 1.0 # days
#fileroot = '_sourcename' 
#meti = 239557417 # value if using "START" for data acquisition
#metf = 255398400 # specific value depending on data downloaded

fermimerge(0)







    


























print('Run time (seconds):')
print (time.time() - start_time)