Python script for ISA atmospheric calculation and velocity conversion

Here is small python script which calculates atmospheric properties and does velocity

conversion: Let me know if you need any explanation.

#!/usr/bin/python2.5

#####################################################################

## Script: Python script to calculate ISA properties and does velocity conversion ##

## ##

## Author: A. Khare ##

## Date: 18.04.2011 ##

## References: ESDU 77022 ##

## ##

#####################################################################

#Import modules

import sys, shutil, os, math

#Define global coefficients

global Gamma, R

Gamma = 1.4

R = 287.05287

#Define atmospheric coefficients

def Coeff(Altitude):

global a, b, c, D, E, F, G, i, j, L, m, n

if Altitude <= 36089.2:

a = 288.15

b = -1.9812 * 10**-3

c = 4.2927085

D = -29.514885 * 10**-6

E = 5.2558797

i = 0.24179285

j = -1.6624675 * 10**-6

L = 4.2558797

elif Altitude >= 36089.2 and Altitude <= 65616.8:

a = 216.65

b = 0

F = 2678.442

G = -48.063462 * 10**-6

m = 4.0012122 * 10**-3

n = -48.063462 * 10**-6

elif Altitude >= 65616.8 and Altitude <= 104987:

a = 196.65

b = 0.3048 * 10**-3

c = 0.79011202

D = 1.2246435 * 10**-6

E = -34.163218

i = 1.1616564

j = 1.8005232 * 10**-6

L = 35.163218

elif Altitude >= 104987 and Altitude <= 154199:

a = 139.05

b = 0.85344 * 10**-3

c = 0.47958369

D = 2.943516 * 10**-6

E = -12.201149

i = 1.3544167

j = 8.3129335 * 10**-6

L = -13.201149

elif Altitude >= 154199 and Altitude <= 164042:

a = 270.65

b = 0

F = 873.60472

G = -38.473855 * 10**-6

m = 1.04466 * 10**-3

n = -38.473855 * 10**-6

#

#Calculate temperature ratio

def Temperature_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Coeff(Altitude)

temp_ratio = (a + b*Altitude)/288.15

return temp_ratio

#Calculate pressure ratio

def Pressure_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Coeff(Altitude)

if b < 0 or b > 0:

Pressure_Ratio = (c + D*Altitude)**E

else:

Pressure_Ratio = F*math.exp(G*Altitude)

Pressure_Ratio = (Pressure_Ratio*4.4482)/(0.3048**2)/101325

return Pressure_Ratio

#Calculate density ratio

def Density_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Density_Ratio = Pressure_Ratio(Altitude)/Temperature_Ratio(Altitude)

return Density_Ratio

#Calculate speed of sound

def Sonic_Velocity(Static_Temp_K):

if Static_Temp_K < 0:

print "Static temperature in Kelvin is out of valid range, Its below zero"

sys.exit(0)

Sonic_Velocity = (Gamma*R*Static_Temp_K)**0.5

return Sonic_Velocity

#Calculate dynamic viscosity ratio

def Dynamic_Viscosity_Ratio(Altitude):

Ts = Temperature_Ratio(Altitude)*288.15

Mu = (1.458 * (10**-6)) * (Ts**1.5) / (Ts + 110.4)

Dynamic_Viscosity_Ratio = Mu / (17.894 * (10**-6))

return Dynamic_Viscosity_Ratio

# Convert KCAS to Mach number

def CAS_to_Mach(Altitude, KCAS):

Ps = Pressure_Ratio(Altitude)

CAS_to_Mach = ((((1 + 0.2 * (KCAS / 661.48)**2)**3.5 - 1) * (1 / Ps) + 1)**(1 / 3.5) - 1)**0.5 * 2.2361

return CAS_to_Mach

# Convert KCAS to KTAS

def CAS_to_TAS(Altitude, KCAS):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

CAS_to_TAS = ((((1+0.2*(KCAS/661.48)**2)**3.5-1)*(1/Ps)+1)**(1/3.5)-1)**0.5*2.2361*(Gamma*R*Ts)**0.5/0.51444

return CAS_to_TAS

# Convert KTAS to Mach number

def TAS_to_Mach(Altitude, KTAS):

Ts = Temperature_Ratio(Altitude)*288.15

TAS_to_Mach = KTAS / (((Gamma * R * Ts) ** 0.5) / 0.51444)

return TAS_to_Mach

# Convert KTAS to KCAS

def TAS_to_CAS(Altitude, KTAS):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

Mach = KTAS / ((Gamma * R * Ts) ** 0.5 / 0.51444)

TAS_to_CAS = (((((Mach ** 2 * 0.2 + 1) ** 3.5 - 1) * (Ps) + 1) ** 0.2857 - 1) * 5) ** 0.5 * 661.48

return TAS_to_CAS

# Convert Mach number to KCAS

def Mach_to_CAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Mach_to_CAS = (((((((1 + 0.2 * Mach ** 2) ** 3.5) - 1) * Ps + 1) ** (1 / 3.5)) - 1) * 5 * 661.48 ** 2) ** 0.5

return Mach_to_CAS

# Convert Mach number to KTAS

def Mach_to_TAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

Mach_to_TAS = Mach * (Gamma * R * Ts) ** 0.5 / 0.51444

return Mach_to_TAS

# Convert Mach number to KEAS

def Mach_to_KEAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

tas = Mach * (Gamma * R * Ts) ** 0.5 / 0.51444

Mach_to_KEAS = tas * Density_Ratio(Altitude) ** 0.5

return Mach_to_KEAS

if len(sys.argv) == 3:

func = sys.argv[1]

para1 = float(sys.argv[2])

if func == "Temperature_Ratio":

Out = Temperature_Ratio(para1)

elif func == "Pressure_Ratio":

Out = Pressure_Ratio(para1)

elif func == "Density_Ratio":

Out = Density_Ratio(para1)

elif func == "Sonic_Velocity":

Out = Sonic_Velocity(para1)

elif func == "Dynamic_Viscosity_Ratio":

Out = Dynamic_Viscosity_Ratio(para1)

else:

print "\033[1;31mType of argument for atmospheric property calculation is not valid. Valid arguments are:\033[1;m"

print "\033[1;31mTemperature_Ratio, Pressure_Ratio, Density_Ratio, Sonic_Velocity, Dynamic_Viscosity_Ratio\033[1;m"

sys.exit(0)

if func == "Sonic_Velocity":

print "%s at %g Kelvin is: %.6f" % (func, para1, Out)

else:

print "%s at %g feet is: %.6f" % (func, para1, Out)

elif len(sys.argv) == 4:

func = sys.argv[1]

para1 = float(sys.argv[2])

para2 = float(sys.argv[3])

if func == "CAS_to_Mach":

Out1 = CAS_to_Mach(para1, para2)

Out2 = "KCAS"

Out3 = "Mach"

elif func == "CAS_to_TAS":

Out1 = CAS_to_TAS(para1, para2)

Out2 = "KCAS"

Out3 = "KTAS"

elif func == "TAS_to_Mach":

Out1 = TAS_to_Mach(para1, para2)

Out2 = "KTAS"

Out3 = "Mach"

elif func == "TAS_to_CAS":

Out1 = TAS_to_CAS(para1, para2)

Out2 = "KTAS"

Out3 = "KCAS"

elif func == "Mach_to_CAS":

Out1 = Mach_to_CAS(para1, para2)

Out2 = "Mach"

Out3 = "KCAS"

elif func == "Mach_to_TAS":

Out1 = Mach_to_TAS(para1, para2)

Out2 = "Mach"

Out3 = "KTAS"

elif func == "Mach_to_EAS":

Out1 = Mach_to_KEAS(para1, para2)

Out2 = "Mach"

Out3 = "KEAS"

else:

print "\033[1;31mType of argument for aerodynamic velocity calculation is not valid. Valid arguments are:\033[1;m"

print "\033[1;31mCAS_to_Mach, CAS_to_TAS, TAS_to_Mach, TAS_to_CAS, Mach_to_CAS, Mach_to_TAS, Mach_to_EAS\033[1;m"

sys.exit(0)

print "%.6f %s at %g feet is: %.6f %s" % (para2, Out2, para1, Out1, Out3)

else:

print "\033[1;31mWrong number of arguments to the script\033[1;m"

print "Use for Atmospheric Properties:"

print "cfd_atm Property_Type Altitude[feets]"

print "Example:"

print "cfd_atm Pressure_Ratio 30000\n"

print "Use for Aerodynamic Velocities:"

print "cfd_atm Converson_Type Altitude[feets] Input[Velocity/Mach]"

print "Example:"

print "cfd_atm CAS_to_TAS 30000 250"

sys.exit(0)

#End of script

#####################################################################

Script can be downloaded from here:

atm_velo_calculator 

Python script to process Cl values in list of directories

Following is the small python script which post process all the directories inside a given directories and store values of Alpha vs CL in a file “Alpha_CL.dat”.

This is helpful when your are running some polar calculations in CFD++.

Use: script common_name_of_directories

Example: If you have AoA polar calculations with directories name as AoA00, AoA02, AoA04 and so on run the following script as:

script AoA

#!/usr/bin/python2.5

#####################################################################

## Script: Python script to post process CL values for all the ##

## cases in given directory ##

## ##

## Author: A. Khare ##

## Date: 03.03.2011 ##

## ##

#####################################################################

#Import modules

import sys, os, glob, commands

#Check for argument. If argument is not given print the use and exit

if len(sys.argv) != 2:

print "Use: get_cl prefix_for_cases"

exit(1)

wdir = os.getcwd()

#Open file for writing CL values

f = open("Alpha_CL.dat", 'w')

#Reading the prefix and finding out all the directories

list = str(sys.argv[1] + "*")

dir_list = glob.glob(list)

#Loop over all the directories and process CL

for i in dir_list:

j = wdir + '/' + i

os.chdir(j)

alpha = float(commands.getoutput("cat infout1f.inp | grep alpha | head -1 | awk \'{print $2}\'"))

cl = float(commands.getoutput("infout1f 1 | tail -21 | head -1 | awk \'{print $9}\'"))

f.write('%.2f\t%.3f\n' % (alpha, cl))

os.chdir(wdir)

f.close()

#End of script

########################################################################