added scripts

.gitignore

@@ -14,4 +14,6 @@
 !/plots/*.asy
 !/plots/*.py
 
-/tmp
\ No newline at end of file
+/tmp
+/x*
+/*.dat
\ No newline at end of file

ring-result

+#!/usr/bin/python
+
+"""
+this scripts analyses output from TLUSTY and collects all important quantities
+output: total flux, eff. temperature, total column mass, disk height
+
+OUTPUT:
+  results are printed to stdout
+
+EXIT CODE:
+  exit code is 0 for converged rings (no errors found); exit code is 1 othervise
+"""
+
+import os
+import sys
+import math
+import numpy as np
+import scipy
+from os import listdir
+from scipy import integrate
+from scipy.integrate import simps
+from collections import OrderedDict
+
+
+def xformat(value):
+    if isinstance(value, (int, long)):
+        return "%d" % value
+    elif isinstance(value, float):
+        return "%E" % value
+    else:
+        return str(value)
+# end of xformat
+
+
+
+
+
+
+script_name = os.path.abspath(__file__)
+script_path = os.path.dirname(script_name)
+
+result = OrderedDict({'converged': 0})
+
+
+
+
+work_dir = sys.argv[1] if (len(sys.argv) >= 2) else script_path
+
+
+# tests
+if (work_dir is None) or not os.path.isdir(work_dir):
+    print "no target directory ({})".format(work_dir)
+    exit()
+
+
+    
+try:
+    # read module.5 to get input parameters
+    mod5 = open(work_dir+'/fort.5', 'r').readlines();
+    mod5_params = mod5[0].strip().split();
+    mod5_modes  = mod5[1].strip().split();
+
+    result['target_temp']    = float(mod5_params[1]);
+    result['target_flux']    = 5.670400e-05*pow(result['target_temp'],4)
+    result['target_sigma']   = float(mod5_params[3]);
+    result['target_gravity'] = float(mod5_params[2]);
+    result['target_height']  = float(open(work_dir+'/../param-H','r').read());
+    
+
+
+    # module.13 read specific intensities and calculate resulting effective temperature
+    if os.path.isfile(work_dir+'/fort.15'): mod13 = open(work_dir+'/fort.13', 'r');
+    arr_f = []; arr_i = []
+    for line in (x for x in mod13 if not x.startswith('#')):
+        values = line.strip().replace("D","E").split()
+        arr_f.append(float(values[0]))
+        arr_i.append(float(values[1]))
+    if (arr_f[0] > arr_f[1]): arr_f.reverse(); arr_i.reverse()  # reverse order of values if necessary
+    mod13_flux = 4*np.pi*integrate.simps(arr_i,arr_f)
+    mod13.close()
+
+    # from module.15 read specific intensities and calculate resulting effective temperature
+    if os.path.isfile(work_dir+'/fort.15'): mod15 = open(work_dir+'/fort.15', 'r');
+    arr_f = []; arr_i = []
+    for line in (x for x in mod15 if not x.startswith('#')):
+        if 'NaN' in line: raise Exception("mod15 contains NaNs")
+        values = line.strip().split()
+        arr_f.append(float(values[0]))
+        arr_i.append(float(values[1]))
+    if (arr_f[0] > arr_f[1]): arr_f.reverse(); arr_i.reverse()  # reverse order of values if necessary
+    mod15_flux = integrate.simps(arr_i,arr_f)
+    mod15.close()
+    
+    # from module.98 read total column density
+    if os.path.isfile(work_dir+'/fort.98'): mod98 = open(work_dir+'/fort.98', 'r').readlines();
+    depth0 = mod98[ 0].strip().replace("D","E").split()
+    depthN = mod98[-2].strip().replace("D","E").split()
+    mod98_colmass = float(depthN[1])
+    mod98_zheight = float(depth0[2])/1.47703573e+6 # scale to grav.radius units (1.47703573e+6 = GM/c^2 for M=10 in cm)
+    
+    result['flux']  = mod15_flux  # integrate Iv to obtain total flux
+    result['flux_mod13']  = mod13_flux  # integrate Iv to obtain total flux
+    result['flux_mod15']  = mod15_flux  # integrate Iv to obtain total flux
+    result['flux_check'] = int((abs(mod15_flux-mod13_flux) / mod13_flux) < 1e-2)
+    result['temp']  = math.pow(result['flux']/5.670400e-05, 0.25);  # from flux get temperature
+    result['sigma'] = mod98_colmass
+    result['height'] = mod98_zheight
+    result['delta_temp'] = abs(result['temp']-result['target_temp']) / result['target_temp'] 
+    result['delta_sigma'] = abs(result['sigma']-result['target_sigma']) / result['target_sigma'] 
+    result['delta_height'] = abs(result['height']-result['target_height']) / result['target_height'] 
+    
+    result['converged'] = 1
+except:
+    error = '!! Error reading tlusty results in '+work_dir+': ' + str(sys.exc_info()[1])
+    sys.stderr.write(error+'\n')
+    print('# '+error)
+    result['converged'] = 0
+
+
+# print out content of result in key=value format
+for key, value in result.items(): print key+'='+xformat(value)
+
+
+
+sys.exit( 0 if (result['converged']==1) else 1 )
+
+

ring-st1-calc-lte

+#!/bin/bash
+
+# --- set working directory either
+# choose either dir given as argument or the script's home directory
+DIR=$(dirname "$0")
+if [ "$1" != "" ]; then DIR=$1; fi
+
+# --- test if directory exists and if it contains TLUSTY executable
+if [ ! -e $DIR/t205 ]; then
+    echo "ERROR: working dir ($DIR) does not exist or it does not have TLUSTY"
+    exit 1
+fi
+
+# change cwd to $DIR
+cd $DIR
+
+
+# --- test if there is anything to be done (ring is converged already)
+./ring-result >/dev/null 2>/dev/null
+if [ $? -eq 0 ]; then
+    echo "ring OK ... skipping"
+    exit 0
+fi
+
+
+
+# --- setup parameters
+rm -f fort.* 
+TEMP=$(cat ../param-T)
+GRAV=$(cat ../param-Q)
+SIGM=$(cat ../param-S)
+MODEL_LTE="T"
+MODEL_LTEGRAY="T"
+FREQ_POINTS=200
+MAXITER=1
+
+
+# --- setup init file (UNIT 1 and 5)
+rm -f stage1-init
+echo 1 > fort.1
+echo "0.0000  $TEMP  $GRAV  $SIGM" >> stage1-init
+echo "$MODEL_LTE  $MODEL_LTEGRAY" >> stage1-init
+echo "'stage1-options'   'fort'" >> stage1-init
+echo "$FREQ_POINTS" >> stage1-init
+cat t205-atomicdata >> stage1-init
+echo "0.0000  $TEMP  $GRAV  $SIGM" >> stage1-init
+cp stage1-init fort.5
+
+
+
+# --- setup init options file
+rm -f stage1-options
+echo "ALPHAV=0.100,NFTAIL=-41,DFTAIL=0.15,NITZER=11" >> stage1-options
+echo "ICOMPT=7,ICOMST=0,IACC=14,ITEK=40,IPFREC=0" >> stage1-options
+
+
+
+# --- run TLUSTY
+# first, run with gray model anzatz
+nice -n 15 ./t205 < stage1-init > fort.6 2>fort.err
+
+# second, run again providing first run result as an input gray model
+sed -i 's/T  T/T  F/g' stage1-init
+cp stage1-init fort.5
+cp fort.7 fort.8
+nice -n 15 ./t205 < stage1-init > fort.6 2>fort.err
+
+
+
+./ring-result > results
+if [ $? -ne 0 ]; then
+    echo "NO CONVERGENCE"
+    exit 1
+fi
+
+
+
+exit 0
+
+

setup-stage1

+#!/usr/bin/perl
+use File::Basename;
+use File::Path;
+use File::Copy;
+use Cwd 'abs_path';
+
+if ($#ARGV != 1) { die("setup: spin luminosity"); }
+
+my $spin  = $ARGV[0];
+my $lum   = $ARGV[1];
+my $alpha = 0.100;
+my $rmax  = 2000.0;
+my $nrings = 100;
+
+print "# spin    = $spin\n";
+print "# lum     = $lum\n";
+print "# alpha   = $alpha\n";
+print "# rmax    = $rmax\n";
+print "# N-rings = $nrings\n";
+
+my $base_dir = sprintf("a%03.0f/l%03.0f", $spin*1000, $lum*100);
+my $tlusty_dir = 'tlusty';                 # folder with tlusty 
+my $radial_dir = 'disk-model-sd';          # folder with radial model
+
+if (-e $base_dir) { die("Base directory exists, choose a new one!"); }
+if (!(-e $tlusty_dir)) { die("Tluty directory does not exists!"); }
+
+# make $base_dir folder and model index file
+mkpath($base_dir);
+system("(cd $radial_dir; ./radial-model $spin $lum $alpha $rmax $nrings) > $base_dir/index.dat");
+
+# open index file for reading
+open(my $fh, '<', "$base_dir/index.dat")
+    or die("cannot open < radial model : $!");
+
+
+while (<$fh>) {
+    ($N, $R, $T, $S, $Q, $H, $L, $V) = split(/\s+/);
+
+    # make folder     
+    $path = sprintf("%s/%04d", $base_dir, $N);
+    mkpath("$path");
+    mkpath("$path/stage1");
+
+    # copy setup file in
+    symlink(abs_path("$tlusty_dir/t205"), "$path/stage1/t205");
+    symlink(abs_path("$tlusty_dir/t205-atomicdata"), "$path/stage1/t205-atomicdata");
+    symlink(abs_path("$tlusty_dir/t205-options"), "$path/stage1/t205-options");
+    symlink(abs_path("$tlusty_dir/data"), "$path/stage1/data");
+    for $f (glob "$tlusty_dir/*.FOR") { copy($f, "$path/stage1") };
+
+    system("echo $N > $path/param-id");
+    system("echo $R > $path/param-R");
+    system("echo $T > $path/param-T");
+    system("echo $S > $path/param-S");
+    system("echo $Q > $path/param-Q");
+    system("echo $H > $path/param-H");
+    
+    system("rm -f $path/params");
+    system("echo 'ring-id=$N' >> $path/params");
+    system("echo 'radius=$R' >> $path/params");
+    system("echo 'temp=$T' >> $path/params");
+    system("echo 'sigma=$S' >> $path/params");
+    system("echo 'gravity=$Q' >> $path/params");
+    system("echo 'height=$H' >> $path/params");
+    system("echo 'angmom=$L' >> $path/params");
+    system("echo 'radvel=$V' >> $path/params");
+
+    symlink(abs_path("./ring-result"), "$path/stage1/ring-result");
+    symlink(abs_path("./ring-st1-calc-lte"), "$path/stage1/ring-st1-calc-lte");
+#    #symlink(abs_path("./ring-st2-lteg"), "$path/ring-stage2-lteg");
+}
+
+close($fh);
+
+

stage1-check

+#!/usr/bin/perl
+use File::Basename;
+use File::Path;
+use Cwd 'abs_path';
+use File::Slurp;
+use List::Util qw[min max];
+
+if ($#ARGV != 0) { die("stage1-check: dir"); }
+
+my $base_dir  = $ARGV[0];
+
+my @dirs = grep { -d } glob "$base_dir/????";
+
+
+open(output, ">", "stage1-check.dat") or die $!;
+print(output "# checking $base_dir\n");
+print(output "# format:  id | radius | T(sd) | T(tl) | H(sd) | H(tl)\n");
+print(output "# gnuplot: plot 'stage1-check.dat' u 1:(\$3/1e6) w l, '' u 1:(\$4/1e6) w p lw 2 lc 3, '' u 1:5 w l, '' u 1:6 w l\n");
+
+foreach $dir (@dirs) {
+    my $id = trim(read_file("$dir/param-id"));
+    my $R  = trim(read_file("$dir/param-R")); 
+    my $T  = trim(read_file("$dir/param-T")); 
+    my $H  = trim(read_file("$dir/param-H")); 
+
+    my @mod15 = read_file("$dir/stage1/fort.15");
+    my $tlH = trim($mod15[11])/1.476716e+06;
+    my $tlF = 0.0;
+    foreach my $n (12..$#mod15-1) {
+        (my $f0, my $Iv0) = split(/\s+/, trim($mod15[$n]));
+        (my $f1, my $Iv1) = split(/\s+/, trim($mod15[$n+1]));
+        $tlF += exp(0.5*(log($Iv0+1e-200)+log($Iv1+1e-200)) + 0.5*(log($f1)+log($f0))) * abs(log($f1)-log($f0))
+    }
+    my $tlT = sqrt(sqrt($tlF/5.670400e-05));
+
+    #my @mod9 = read_file("$dir/stage1/fort.9");
+    #my $n = $#mod9-1;
+    #my $last_iter = $1 if ($mod9[$n] =~ m/^\s+(\d+)\s/);
+    #my $max_change = 1e-99;
+    #while ($mod9[$n] =~ m/\s+$last_iter\s+.*?\s+.*?\s+.*?\s+.*?\s+(.*?)\s+/i) {
+    #    my $tmp = $1;
+    #    $tmp =~ s/d/e/i;
+    #    $max_change = max($max_change, abs($tmp));
+    #    $n--;
+    #}
+    #print "$dir - $max_change\n";
+
+    printf(output "%04d  %6.2f  %.4e  %.4e  %.4e  %.4e\n", $id,  $R,  $T,  $tlT,  $H,  $tlH);
+}
+close(output);
+
+
+sub trim {
+    return $_[0] =~ s/^\s+|\s+$//rg;
+}
+
+

stage1-invalidate

+#!/bin/bash
+
+echo "NaN" > $1/stage1.15

stage1-run-ir

+#!/bin/bash
+
+DIR=$1
+BDIR=$(pwd)
+
+echo $DIR
+
+for i in `ls -d $DIR/0*`; do
+    cd $BDIR/$i
+    echo "$i - $(pwd)"
+    ./run-lte-ir &
+    cd $BDIR
+    
+    sleep 1
+done
\ No newline at end of file

stage1-run-lte0

+#!/bin/bash
+
+DIR=$1
+if [ "$DIR" == "" ]; then
+    echo "Usage: $0 <directory>"
+    exit 1
+fi
+
+if [ ! -e $DIR/index.dat ]; then
+    echo "Cannot find index.dat in directory $DIR"
+    exit 1
+fi
+
+
+echo "Running: $DIR:"
+
+for i in `ls -d $DIR/????`; do
+    echo $i
+    $i/stage1/ring-st1-calc-lte >/dev/null &
+    sleep 1
+done
\ No newline at end of file

stage1-run-lteg

+#!/bin/bash
+
+DIR=$1
+echo "Running: $DIR:"
+
+for i in `ls -d $DIR/????`; do
+    $i/ring-stage1-lteg &
+    #sleep 1
+done
\ No newline at end of file