Merge pull request #313 from agnwinds/import

Capability to Import Models

.travis.yml

@@ -56,3 +56,4 @@ script:
     - py -i ulx1
     - py -i ngc5548
     - py -i lamp_post
+    - py -i cv_import

README.md

@@ -1,10 +1,12 @@
 # Python
 
 Python is a (confusingly named) Monte Carlo radiative transfer code which uses the Sobolev approximation.
-It has been developed by Knox Long, Christian Knigge, Stuart Sim, Nick Higginbottom and James Matthews. 
+It has been developed by Knox Long, Christian Knigge, Stuart Sim, Nick Higginbottom, James Matthews and Sam Mangham. 
+The code has been used for a variety of research projects invovling the winds of cataclysmic variables, of young stellar objects, and of
+AGN.
 
-The code is not yet ready for use and should not be installed yet. If you are interested in using
-Python please contact [email protected].
+The code is under active development, but we are looking for beta users to test the code, and potentially use it 
+for their own research. If you are interested in using Python please contact [email protected].
 
 ## Travis Build Status
 

examples/beta/cv_import.pf

@@ -0,0 +1,63 @@
+System_type(0=star,1=binary,2=agn)   1
+disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk)   1
+Number.of.wind.components         1
+Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell)   10
+Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var)   1
+File.with.model2read    cv_standard.import.txt
+Wind.dim.in.x_or_r.direction       30
+Wind.dim.in.z_or_theta.direction   30
+disk.atmosphere(0=no,1=yes)  0 
+Atomic_data   data/standard78
+photons_per_cycle   1000000
+Ionization_cycles   15
+spectrum_cycles   10
+Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow,8=matrix_bb,9=matrix_pow)   8
+Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering)   5
+Thermal_balance_options(0=everything.on,1=no.adiabatic)   0
+System_type(0=star,1=binary,2=agn)   1
+Star_radiation(y=1)       1
+Disk_radiation(y=1)       1
+Boundary_layer_radiation(y=1)   0
+Wind_radiation(y=1)   	  1
+Rad_type_for_star(0=bb,1=models)_to_make_wind   0
+Rad_type_for_disk(0=bb,1=models)_to_make_wind   0
+mstar(msol)   0.8
+rstar(cm)   7e+08
+tstar   40000
+msec(msol) 0.6
+period(hr) 5.57
+disk.mdot(msol/yr)   1e-8
+Disk.illumination.treatment(0=no.rerad,1=high.albedo,2=thermalized.rerad,3=analytic)   0
+Disk.temperature.profile(0=standard;1=readin)   0
+disk.radmax(cm)   2.4e+10
+wind.radmax(cm)   1e+12
+wind.t.init   40000
+wind.mdot(msol/yr)   1e-9
+sv.diskmin(wd_rad)   4
+sv.diskmax(wd_rad)   12
+sv.thetamin(deg)   20
+sv.thetamax(deg)   65
+sv.mdot_r_exponent   0
+sv.v_infinity(in_units_of_vescape   3
+sv.acceleration_length(cm)   7e10
+sv.acceleration_exponent   1.5
+filling_factor(1=smooth,<1=clumpted)		1
+Rad_type_for_star()_in_final_spectrum   0
+Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum   0
+spectrum_wavemin   800
+spectrum_wavemax   1850
+no_observers   5
+angle(0=pole)   10
+angle(0=pole)   27.5
+angle(0=pole)   45
+angle(0=pole)   62.5
+angle(0=pole)   80
+phase(0=inferior_conjunction)   0.5
+phase(0=inferior_conjunction)   0.5
+phase(0=inferior_conjunction)   0.5
+phase(0=inferior_conjunction)   0.5
+phase(0=inferior_conjunction)   0.5
+live.or.die(0).or.extract(anything_else)   1
+spec.type(flambda(1),fnu(2),basic(other)   1
+reverb.type(0=None,1=Photon,2=Wind) 	   0
+Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic)    2

examples/regress/rtheta.pf

@@ -0,0 +1,68 @@
+System_type(0=star,1=binary,2=agn,3=previous)   1
+mstar(msol)                                0.8
+rstar(cm)                                  7e+08
+Star_radiation(y=1)                        1
+Rad_type_for_star(0=bb,1=models)_to_make_wind   0
+tstar                                      40000
+msec(msol)                                 0.6
+period(hr)                                 5.57
+disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk)   1
+Disk_radiation(y=1)                        1
+Rad_type_for_disk(0=bb,1=models)_to_make_wind   0
+disk.mdot(msol/yr)                         1e-8
+Disk.temperature.profile(0=standard;1=readin,2=analytic)   0
+disk.radmax(cm)                            2.4e+10
+Boundary_layer_radiation(y=1)              0
+Wind_radiation(y=1)                        1
+Number.of.wind.components                  1
+Wind_type(0=SV,1=Star,3=Hydro,4=corona,5=knigge,6=homologous,7=yso,9=shell,11=imported)   0
+Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var)   2
+Wind.dim.in.x_or_r.direction               30
+Wind.dim.in.z_or_theta.direction           30
+photons_per_cycle                          100000
+Ionization_cycles                          10
+spectrum_cycles                            10
+Wind_ionization(0=on.the.spot,1=LTE(tr),2=fixed,3=recalc_bb,4=LTE(t_e),6=pairwise_bb,7=pairwise_pow,8=matrix_bb,9=matrix_pow)   8
+Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering)   5
+Atomic_data                                data/standard78
+Surface.reflection.or.absorption(0=no.rerad,1=high.albedo,2=thermalized.rerad)   0
+Thermal_balance_options(0=everything.on,1=no.adiabatic)   0
+wind.radmax(cm)                            1e+12
+wind.t.init                                40000
+wind.mdot(msol/yr)                         1e-9
+sv.diskmin(units_of_rstar)                 4
+sv.diskmax(units_of_rstar)                 12
+sv.thetamin(deg)                           20
+sv.thetamax(deg)                           65
+sv.mdot_r_exponent                         0
+sv.v_infinity(in_units_of_vescape          3
+sv.acceleration_length(cm)                 7e10
+sv.acceleration_exponent                   1.5
+filling_factor(1=smooth,<1=clumped)        1
+Rad_type_for_star(0=bb,1=models,2=uniform)_in_final_spectrum   0
+Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum   0
+spectrum_wavemin                           800
+spectrum_wavemax                           1850
+no_observers                               5
+angle(0=pole)                              10
+angle(0=pole)                              27.5
+angle(0=pole)                              45
+angle(0=pole)                              62.5
+angle(0=pole)                              80
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+live.or.die(0).or.extract(anything_else)   1
+spec.type(flambda(1),fnu(2),basic(other)   1
+reverb.type                                0
+Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic)   2
+ 0.5
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+phase(0=inferior_conjunction)              0.5
+live.or.die(0).or.extract(anything_else)   1
+spec.type(flambda(1),fnu(2),basic(other)   1
+reverb.type                                0
+Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic)   2

py_progs/import_1d.py

@@ -0,0 +1,160 @@
+#!/usr/bin/env python 
+
+'''
+                    Space Telescope Science Institute
+
+Synopsis:  
+
+Read the master file produced by windwave2table for a
+stellar-type model and produce a file to be read by 
+import in python
+
+
+Command line usage (if any):
+
+    usage: import_1d.py filename
+
+Description:  
+
+Primary routines:
+
+    doit
+
+Notes:
+                                       
+History:
+
+171106 ksl Coding begun
+
+'''
+
+import sys
+from astropy.io import ascii
+import numpy
+
+
+def read_file(filename,char=''):
+    '''
+    Read a file and split it into words, eliminating comments
+    
+    char is an optional parameter used as the delimiter for
+    splitting lines into words.  Otherwise white space is
+    assumed.
+
+    History:
+    
+    110729    ksl    Added optional delimiters
+    141209    ksl    Reinstalled in my standard startup
+            script so there was flexibility to
+            read any ascii file
+    '''
+
+    try:
+        f=open(filename,'r')
+        xlines=f.readlines()
+        f.close()
+    except IOError :
+        print ("The file %s does not exist" % filename)
+        return []   
+
+    lines=[]
+
+    i=0
+    while i<len(xlines):
+        z=xlines[i].strip()
+        if char=='':
+            z=z.split()
+        else:
+            z=z.split(char)
+        if len(z)>0:
+            if z[0][0]!='#':
+                lines=lines+[z]
+        i=i+1
+    return lines
+
+
+
+
+def read_table(filename='foo.txt',format=''):
+    '''
+    Read a file using astropy.io.ascii and 
+    return this 
+
+    Description:
+
+    Notes:
+
+    History:
+
+
+    '''
+    try:
+        if format=='':
+            data=ascii.read(filename)
+        else:
+            data=ascii.read(filename,format=format)
+        for col in data.itercols():
+            if col.dtype.kind in 'SU':
+                data.replace_column(col.name,col.astype('object'))
+    except IOError:
+        print ('Error: file %s does not appear to exist' % filename)
+        return
+
+    print ('Here are the column names:')
+
+    print (data.colnames)
+
+    return data
+
+
+def doit(root='star',outputfile=''):
+    '''
+    Do something magnificent
+
+    Description:
+
+    Notes:
+
+    History:
+
+    '''
+
+    filename=root+'.0.master.txt'
+    if outputfile=='':
+        outputfile=root+'.import.txt'
+
+
+
+    data=read_table(filename)
+
+    v=data['v_x']**2+data['v_y']**2+data['v_z']**2
+    v=numpy.sqrt(v)
+
+    data['v']=v
+
+    xdata=data['i','r','v','rho','t_e']
+
+
+    print (xdata)
+
+
+    # This format is the easy to read back automatically
+    ascii.write(xdata,outputfile,format='fixed_width_two_line')
+
+    return
+
+
+
+
+
+
+
+
+# Next lines permit one to run the routine from the command line
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv)>1:
+        # doit(int(sys.argv[1]))
+        doit(sys.argv[1])
+    else:
+        print ('usage: import_1d.py filename')