Merge branch 'master' of github.com:othocaes/cloudy-agn

scripts/meta/bulk_interpolation_fix.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+
+script_dir=$( cd $(dirname $0) ; pwd -P |sed 's@^\(.*\)/scripts.*@\1/scripts@')
+basedir=`pwd`
+
+if [[ -x $script_dir/../bin/interpolation_fix ]]
+then
+    echo "Binary available. Proceeding."
+else
+    echo "Binary unavailable. Exiting."
+    exit 9
+fi
+
+for fortfile in $@
+do
+    if [[ -e $(basename $fortfile) ]]
+    then
+        mkdir -p uninterpolated
+        mv $fortfile uninterpolated
+        $script_dir/../bin/interpolation_fix uninterpolated/$(basename $fortfile) $fortfile
+    else
+        $script_dir/../bin/interpolation_fix $fortfile $(basename $fortfile)
+    fi
+    echo ""
+done
+exit 0

scripts/meta/create_fort_files_and_interpolate.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+
+# This script not generalized. Use with caution.
+
+mkdir fortfiles
+cd fortfiles
+../../../../../bin/create_fort_files ../mpi_grid.out ../../../../../reference/linelist
+../../../../../scripts/meta/bulk_interpolation_fix.sh fort.* >> interpolation_report
+tar cf ../fortfiles_4thdex_ngc5548_mehdipour_cldn24.tar fort.*
+cd ..
+gzip fortfiles_4thdex_ngc5548_mehdipour_cldn24.tar
+echo done

src/agn.hpp

@@ -45,7 +45,7 @@ struct cloudy_line_data {
 
 typedef std::map<std::string,cloudy_line_data> cloudy_line_output;
 
-struct cloudy_result {	
+struct cloudy_result {
 	std::string header, footer;
 	std::list<std::string> emergent_line_raw_text;
 	std::list<std::string> cautions;
@@ -60,7 +60,7 @@ struct cloudy_result {
 		iterations(0),
 		phi(0),
 		hden(0),
-		colden(0) 
+		colden(0)
 		{}
 };
 
@@ -70,7 +70,7 @@ cloudy_grid read_cloudy_grid(std::ifstream&);
 //	Operator<< prints general info about the run result.
 std::ostream& operator<< (std::ostream&, cloudy_result);
 //  Seeks an instream to the line after the first occurrence of seek_string.
-void seek_to(std::string,std::istream&); 
+void seek_to(std::string,std::istream&);
 
 
 } // end namespace agn
@@ -123,7 +123,7 @@ void agn::seek_to(std::string seek_string,std::istream& stream) {
 
 std::string agn::format_table1d(agn::table1d table) {
     std::stringstream output;
-    output 
+    output
         << std::scientific
         << std::setprecision(5);
     agn::table1d::iterator x=table.begin();
@@ -272,7 +272,7 @@ agn::cloudy_grid agn::read_cloudy_grid(std::ifstream& inputfile) {
 			}
 		}
 		point.colden = colden;
-		if(agn::debug) std::cout 
+		if(agn::debug) std::cout
 			<< " Grabbing emission lines.";
 		agn::seek_to("Emergent line intensities",inputfile);
 		std::list<std::string> emergent_line_raw_text;
@@ -321,7 +321,7 @@ agn::cloudy_grid agn::read_cloudy_grid(std::ifstream& inputfile) {
 			}
 			linetext_it++;
 		}
-		if(agn::debug) std::cout 
+		if(agn::debug) std::cout
 			<< " Grabbing footer.";
 		while (inputline == "")
 			getline(inputfile,inputline);
@@ -332,9 +332,9 @@ agn::cloudy_grid agn::read_cloudy_grid(std::ifstream& inputfile) {
 			footer.append("\n");
 			getline(inputfile,inputline);
 		}
-		if(agn::debug) std::cout 
+		if(agn::debug) std::cout
 			<< "\nAdding point to grid: "
-			<< point 
+			<< point
 			<< std::endl;
 		grid[*coords] = point;
 		coords++;

src/create_fort_files.cpp

@@ -19,8 +19,8 @@ int main(int argc, char const *argv[]) {
 	std::ifstream line_list_file;
 	line_list_file.open(argv[2]);
 	agn::line_list lines_to_print = agn::read_line_list(line_list_file);
-	std::cout 
-		<< "Compiling table2ds for " 
+	std::cout
+		<< "Compiling table2ds for "
 		<< lines_to_print.size()
 		<< " emission lines.\n";
 	std::list<agn::eqwidth_table> tables = agn::compile_eqwidth_tables(grid,lines_to_print,1215.00);
@@ -43,10 +43,46 @@ int main(int argc, char const *argv[]) {
 		fortfilenum++;
 	}
 
+    int num_unconverged = 0;
+    agn::cloudy_grid::iterator result_it = grid.begin();
+    std::ofstream cautionreportfile;
+    cautionreportfile.open("cautions");
+    cautionreportfile
+        << "The following solutions probably did not converge."
+        << std::endl << std::endl;
+    while(result_it != grid.end()) {
+        if (result_it->second.iterations >= 40) {
+            num_unconverged++;
+            cautionreportfile
+                << "hden = "
+                << std::fixed
+                << std::setprecision(3)
+                << result_it->second.hden
+                << ", phi = "
+                << result_it->second.phi
+                << std::endl
+                << "───────────────────────────"
+                << std::endl;
+            std::list<std::string>::iterator caution_it = result_it->second.cautions.begin();
+            while(caution_it != result_it->second.cautions.end()) {
+                cautionreportfile
+                    << *caution_it
+                    << std::endl;
+                    caution_it++;
+            }
+            cautionreportfile << std::endl << std::endl;
+        }
+        result_it++;
+    }
+
+    std::cout
+        << "Saved cautions for "
+        << num_unconverged
+        << " unconverged solutions."
+        << std::endl;
+        cautionreportfile.close();
+
 	std::cout << "Done.\n";
 	return 0;
 }
 
-
-
-

src/interpolation_fix.cpp

@@ -34,7 +34,7 @@ int main(int argc, char const *argv[]) {
             << outlier_it->second
             << "\n"
             << std::scientific
-            << "Old value: " 
+            << "Old value: "
             << jagged_table.value[outlier_it->first][outlier_it->second]
             << "\t"
             << "New value: "
@@ -43,7 +43,7 @@ int main(int argc, char const *argv[]) {
             outlier_it++;
     }
     std::ofstream smoothed_table_file;
-    smoothed_table_file.open("smoothed_eqwidth_table");
+    smoothed_table_file.open(argv[2]);
     smoothed_table_file << agn::format_eqwidth_table(smooth_table);
     smoothed_table_file.close();
 
@@ -118,13 +118,13 @@ agn::eqwidth_table agn::smooth(agn::eqwidth_table & jagged_table,agn::gridcoordl
         if ( agn::debug ) {
             std::cout
                 << std::scientific
-                << "Smoothed: " 
+                << "Smoothed: "
                 << newx
                 << ", "
-                << oldy 
+                << oldy
                 << " -> "
                 << newy
-                << "\n"; 
+                << "\n";
         }
 
         coord_it++;
@@ -150,23 +150,23 @@ agn::gridcoordlist agn::known_outliers() {
 }
 
 agn::gridcoordlist agn::find_outliers(agn::eqwidth_table jagged_table) {
-    if (agn::debug) std::cout 
+    if (agn::debug) std::cout
         << "Scanning for outliers.\n";
     agn::gridcoordlist outliers;
 
-    // The distribution is 
-    // parameterised by x=log(jagged_table.hden) and 
+    // The distribution is
+    // parameterised by x=log(jagged_table.hden) and
     // y=log(jagged_table.eqwidth).
 
     agn::iterator2d hden_it = jagged_table.value.begin();
     while ( hden_it != jagged_table.value.end() ) {
         agn::table1d slice = hden_it->second;
-        if( agn::debug ) std::cout 
+        if( agn::debug ) std::cout
             << "\nhden= "
-            << (hden_it->first) 
+            << (hden_it->first)
             << ": ";
 
-        // Crawl, checking slope, until a positive slope is found. 
+        // Crawl, checking slope, until a positive slope is found.
         agn::table1d::iterator phi_it = slice.begin();
         double x1,x2,y1,y2,slope;
         std::vector<double> ref_curve_x,ref_curve_y,curve_back_x,curve_back_y;
@@ -220,7 +220,7 @@ agn::gridcoordlist agn::find_outliers(agn::eqwidth_table jagged_table) {
             phi_it++;
         }
 
-        if ( phi_it == slice.end() ) { 
+        if ( phi_it == slice.end() ) {
             if (agn::debug) std::cout
                 << "Clean";
             hden_it++;
@@ -282,7 +282,7 @@ agn::gridcoordlist agn::find_outliers(agn::eqwidth_table jagged_table) {
 
         std::vector<double>::iterator curve_back_x_it = curve_back_x.begin();
         std::vector<double>::iterator curve_back_y_it = curve_back_y.begin();
-        while ( curve_back_x_it != curve_back_x.end() && 
+        while ( curve_back_x_it != curve_back_x.end() &&
                  curve_back_y_it != curve_back_y.end() ) {
             ref_curve_x.push_back(*curve_back_x_it);
             ref_curve_y.push_back(*curve_back_y_it);
@@ -329,7 +329,7 @@ agn::gridcoordlist agn::find_outliers(agn::eqwidth_table jagged_table) {
         //    << ", ";
         phi_it = slice.find(anomalies_start_x);
         while (phi_it != (++(slice.find(anomalies_end_x)))) {
-            y1 = log10(phi_it->second); 
+            y1 = log10(phi_it->second);
             x1 = phi_it->first;
             y2 = testing_spline[x1];
             //threshold = RATIO_THRESHOLD_MULTIPLIER * (y2 - EQWIDTH_MIN_VAL_LOG);

src/spectral_lines.hpp

@@ -143,7 +143,7 @@ std::string agn::format_eqwidth_table_slice(agn::eqwidth_table table,agn::iterat
 
 
 std::ostream& agn::operator<< (std::ostream& outstream, agn::eqwidth_table table) {
-	outstream 
+	outstream
 		<< table.header[0]
 		<< std::endl
 		<< table.header[1]
@@ -211,9 +211,9 @@ std::list<agn::eqwidth_table> agn::compile_eqwidth_tables(agn::cloudy_grid grid,
 		double x,y;
 		agn::cloudy_line_data data;
 		agn::eqwidth_table new_table;
-        if(agn::debug) std::cout 
-            << "Processing label " 
-            << *line_label_it 
+        if(agn::debug) std::cout
+            << "Processing label "
+            << *line_label_it
             << std::endl;
 		new_table.header[0] = *line_label_it;
 		new_table.header[0].append("  relative to Inci 1215.00A scaled to 1215.00A");