Saving ramses outputs in gasoline format to use with AHF

ramses_to_AHF.py

def prepare_for_amiga(outname, write = False, run_amiga=False) :
    """
    Takes a RAMSES output and turns it into a 'tipsy' file for use with 
    tipsy, pkdgrav, or AHF
    """
    
    import os 
    from pynbody.units import Unit

    s = pynbody.load(outname)
        
    # figure out the units starting with mass

    cmtokpc = 3.2407793e-22
    lenunit  = s._info['unit_l']/s.properties['a']*cmtokpc
    massunit = pynbody.analysis.cosmology.rho_crit(s,z=0,unit='Msol kpc^-3')*lenunit**3
    G_u = 4.4998712e-6 # G in kpc^3 / Msol / Gyr^2
    timeunit = np.sqrt(1/G_u * lenunit**3/massunit)
    
    l_unit = Unit('%f kpc'%lenunit)
    t_unit = Unit('%f Gyr'%timeunit)
    v_unit = l_unit/t_unit
   
    print massunit, timeunit

    newfile = "%s_tipsy/%s_fullbox.tipsy"%(s.filename,outname)

    if write:
        s['mass'].convert_units('%f Msol'%massunit)
        s.g['temp']
        print s['mass']
        s.g['metals'] = s.g['metal']
        s['pos'].convert_units(l_unit)
        s['vel'].convert_units(v_unit)
        s['eps'] = s.g['smooth'].min()
        s['eps'].units = s['pos'].units
        del(s.g['metal'])
        del(s['smooth'])
        print s['vel']
        print s['pos']
        
        s.write(filename='%s'%newfile, fmt=pynbody.tipsy.TipsySnap, binary_aux_arrays = True)

    if run_amiga : spawn_amiga(s,newfile,lenunit, massunit, timeunit)

def spawn_amiga(s, newfile, lenunit, massunit, timeunit, zbox = False) :
    """
    Writes an AHF parameter file with the (hopefully) correct system of units
    and runs the halo finder.
    """
    
    from pynbody.units import Unit, G
    import os
    
    l_unit = Unit('%f kpc'%lenunit)
    t_unit = Unit('%f Gyr'%timeunit)
    v_unit = l_unit/t_unit
   
    f = open('%s.AHF.input'%newfile,'w')
    f.write('[AHF]\n')
    f.write('ic_filename = %s\n'%newfile)
    f.write('ic_filetype = 90\n')
    f.write('outfile_prefix = %s\n'%newfile)
    f.write('LgridDomain = 256\n')
    f.write('LgridMax = 2097152\n')
    f.write('NperDomCell = 5\n')
    f.write('NperRefCell = 5\n')
    f.write('VescTune = 1.0\n')
    f.write('NminPerHalo = 50\n')
    f.write('RhoVir = 0\n')
    f.write('Dvir = 200\n')
    f.write('MaxGatherRad = 1.0\n')
    f.write('[TIPSY]\n')
    f.write('TIPSY_BOXSIZE = %e\n'%(s.properties['boxsize'].in_units('Mpc')*s.properties['h']/s.properties['a']))
    f.write('TIPSY_MUNIT   = %e\n'%(massunit*s.properties['h']))
    f.write('TIPSY_OMEGA0  = %f\n'%s.properties['omegaM0'])
    f.write('TIPSY_LAMBDA0 = %f\n'%s.properties['omegaL0'])
    
 #   velunit = Unit('%f cm'%s._info['unit_l'])/Unit('%f s'%s._info['unit_t'])
    
    f.write('TIPSY_VUNIT   = %e\n'%v_unit.ratio('km s^-1 a', **s.conversion_context()))
    

    # the thermal energy in K -> km^2/s^2

    f.write('TIPSY_EUNIT   = %e\n'%((pynbody.units.k/pynbody.units.m_p).in_units('km^2 s^-2 K^-1')*5./3.))
    f.close()
    
    else : 
        os.environ['OMP_NUM_THREADS'] = '16'
        os.system("~/bin/amiga_pthread_for_tipsyramses %s.AHF.input"%newfile)