TIMD integration.py

import math, sys
import os.path
from math import sqrt

#to use the script, choose the right filename prefix
#################################################################################################################
## loops over "fn" of output files .
## fn = file number
## lambda is the name of a function in python!!!!don't use it as a variable. 
 
def read_dvdl_out(filename):

    lamb=[];dvdl = [];err = [];
    file_op = open(filename, 'r');

    lines  =  file_op.readlines();

    for line in lines:
        linesplit = line.split() #split on white space
        if ( len(linesplit) > 0 ):
            lamb.append(float(linesplit[0]));
            dvdl.append(float(linesplit[1]));
            err.append(float(linesplit[2]));
    file_op.close();

    if ( lamb[0] != 0 ):
        ynew = (lamb[0]*dvdl[1]-lamb[1]*dvdl[0])/(lamb[0]-lamb[1]);
        lamb.insert(0, float(0));
        dvdl.insert(0, ynew);
        err.insert(0, err[0]);

    if ( lamb[len(lamb)-1] != 1):
        i = len(lamb);
        ynew = ( dvdl[i-1]*(lamb[i-2]-1)+dvdl[i-2]*(1-lamb[i-1]) )/(lamb[i-2]-lamb[i-1]);
        lamb.append(float(1));
        dvdl.append(ynew);
        err.append(err[i-1]);

    i = len(lamb);
    ##print i;
    width = [];
    for j in range(0,i):
        if (i == 1 ):
            width.insert(j,1);
        elif (j == 0 ):
            width.insert(j,0.5 * (lamb[j] + lamb[j+1]));
        elif (j == i-1):
            width.insert(j,1-0.5*(lamb[j]+lamb[j-1]));
        else:
            width.insert(j,0.5*(lamb[j+1]-lamb[j-1]));

    return lamb,dvdl,err,width;
#################################################################################################################
def integration_dvdl( lamb,dvdl,err,width ):
    tot_dvdl = 0;
    tot_err = 0;
    i = len(lamb);
    for j in range(0,i):
        tot_dvdl += width[j]*dvdl[j];
        tot_err += width[j]*err[j];
    return tot_dvdl, tot_err;
    
######
def write_dvdl_prep( dvdl,rms_dvdl,file):

    length = len(dvdl);
    for i in range(1,length+1):
        file.write( "%5.3f%9.3f%6.3f" %(i/10.0,dvdl[i-1],rms_dvdl[i-1]) );
        file.write( "\n" );
    return;
#################################################################################################################
def main():

########integration for dvdl at different lambda values
    filename = "integ_prep_step1";
    lamb,dvdl,err,width = read_dvdl_out(filename);
    tot_dvdl, tot_err = integration_dvdl(lamb,dvdl,err,width);
    print "Total DV/DL:",tot_dvdl,"+-",tot_err;

#################################################################################################################
main()