20140725

======================================================

#! /usr/bin/python
#
import math
import os
import time
#
# Run a series of Radex models to estimate temperature & density
# from observed ratios of HCO+ 1-0/3-2, 3-2/4-3, and 4-3/7-6 lines.
#
# Floris van der Tak, 27oct05; revised 28jun06
#
# Grid boundaries
#
tmin = 10.0  # minimum kinetic temperature (K)
tmax = 200.0 # maximum kinetic temperature (K)
nmin = 1e3   # minimum H2 density (cm^-3)
nmax = 1e8   # maximum H2 density (cm^-3)
#
# Parameters to keep constant
#
tbg   = 2.73 # background radiation temperature
cdmol = 1e12 # low enough to be optically thin
dv    = 1.0  # line width (km/s)
#
# Numerical parameters
#
ntemp = 30  # number of temperature points
ndens = 30  # number of density points
bw    = 0.001 # "bandwidth": free spectral range around line
#
# No user changes needed below this point.
#
#!!!!!!!!!!!!!!!!!!!!!!!!
def write_input(infile,tkin,nh2):
   infile.write(mole+'.datn')
   infile.write('radex.outn')
   infile.write(str(flow*(1-bw))+' '+str(fupp/(1-bw))+'n')
   infile.write(str(tkin)+'n')
   infile.write('1n')
   infile.write('H2n')
   infile.write(str(nh2)+'n')
   infile.write(str(tbg)+'n')
   infile.write(str(cdmol)+'n')
   infile.write(str(dv)+'n')
#!!!!!!!!!!!!!!!!!!!!!!!!!!!

def read_radex(outfile):
   line  = outfile.readline()
   words = line.split()
   while (words[1] != "T(kin)"):
       line  = outfile.readline()
       words = line.split()
   temp  = float(words[-1])
   line  = outfile.readline()
   words = line.split()
   dens  = float(words[-1])
   while (words[-1] != "FLUX"):
       line  = outfile.readline()
       words = line.split()
   line  = outfile.readline()
   line  = outfile.readline()
   words = line.split()
   ftmp  = float(words[4])
   while ((ftmp < flow*(1-bw)) or (ftmp > flow/(1-bw))):
       line  = outfile.readline()
       words = line.split()
       ftmp  = float(words[4])
   low   = float(words[-2])
   line  = outfile.readline()
   words = line.split()
   ftmp  = float(words[4])
   while ((ftmp < fupp*(1-bw)) or (ftmp > fupp/(1-bw))):
       line  = outfile.readline()
       words = line.split()
       ftmp  = float(words[4])
   upp   = float(words[-2])
   ratio = low / upp
   return temp,dens,ratio
 
# Begin of main program

start = time.time()

mole = 'hco+'
#!!!!!!!!!!!!!!!!!!!!!
# frequency of the energy level and file name
acts = ([89.2,267.6,'1-0_3-2.dat'],[267.6,356.7,'3-2_4-3.dat'],[356.7,624.2,'4-3_7-6.dat'])
#!!!!!!!!!!!!!!!!!!!!!

for act in acts:
   flow = act[0]
   fupp = act[1]
   gfil = act[2]
#!!!!!!!!!!!!!!!!!!!!!
# input file
   infile = open('radex.inp','w')
#!!!!!!!!!!!!!!!!!!!!!
   print "Starting ",gfil

   for itemp in range(ntemp+1):
       for idens in range(ndens+1):

           temp = tmin*((tmax/tmin)**(float(itemp)/ntemp))
           dens = nmin*((nmax/nmin)**(float(idens)/ndens))
           write_input(infile,temp,dens)
           if (itemp == ntemp and idens == ndens):
               infile.write('0n')
               infile.close()
           else:
               infile.write('1n')

   os.system('radex < radex.inp > /dev/null')
   grid = open(gfil,'w')
#!!!!!!!!!!!!!!!!!!!!!!!!!!!
# format
   fmt  = '%10.3e %10.3e %10.3e n'
#!!!!!!!!!!!!!!!!!!!!!!!!!!!
   outfile  = open('radex.out')

   rmin = 100
   rmax = 0.1

   for itemp in range(ntemp+1):
       for idens in range(ndens+1):

           temp = tmin*((tmax/tmin)**(float(itemp)/ntemp))
           dens = nmin*((nmax/nmin)**(float(idens)/ndens))
           temp,dens,ratio = read_radex(outfile)

           if (ratio > 0.0):
               if (ratio < rmin):
                   rmin = ratio
               if (ratio > rmax):
                   rmax = ratio
#!!!!!!!!!!!!!!!!!!!!!!!!!!
# output
           grid.write(fmt %(temp, math.log10(dens), ratio))
#!!!!!!!!!!!!!!!!!!!!!!!!!!
   grid.close()
   outfile.close()

   print "Min, max:", rmin, rmax

stop = time.time()
dure = stop - start
print "Run time = ",dure, "seconds"

======================================================

20140628

Radex的目录bin里提供了两个python脚本
radex_column.py可以根据提过的参数计算线强
radex_grid.py计算一系列温度和压强值下的线强比

20140627

输入参数文件各项的含义如下（其实和网页版Radex是一样的）

hco+.dat   file containing molecular data (see here)
hco+.rdx   file to write output to (will overwrite)
50 500     output frequency range (GHz; 0 0 means unlimited)
20         kinetic temperature (K)
2          number of collision partners (see under Details)
H2         first collision partner
10000      density of first collision partner (cm-3)
e          second collision partner
1          density of second collision partner (cm-3)
2.73       temperature of background radiation (K; see under Details)
1e12       molecular column density (cm-2)
2.0        line width (km/s)
1 (or 0)   run another calculation (or not)

20140626

之前已经写过一点Radex的使用方法了。

http://blog.sciencenet.cn/blog-117333-733025.html

打算根据CS线强比计算一下密度和温度。