00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ParseDLaw parse parameters on the dlaw command */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 #include "optimize.h" 00007 #include "abund.h" 00008 #include "input.h" 00009 #include "radius.h" 00010 #include "parser.h" 00011 00012 void ParseDLaw(Parser &p) 00013 { 00014 bool lgEnd; 00015 long int j; 00016 00017 DEBUG_ENTRY( "ParseDLaw()" ); 00018 00019 if( dense.gas_phase[ipHYDROGEN] > 0. ) 00020 { 00021 fprintf( ioQQQ, " PROBLEM DISASTER More than one density command was entered.\n" ); 00022 cdEXIT(EXIT_FAILURE); 00023 } 00024 00025 /* call fcn dense_fabden(RADIUS) which uses the ten parameters 00026 * N.B.; existing version of dense_fabden must be deleted 00027 * >>chng 96 nov 29, added table option */ 00028 if( p.nMatch("TABL") ) 00029 { 00030 /* when called, read in densities from input stream */ 00031 strcpy( dense.chDenseLaw, "DLW2" ); 00032 if( p.nMatch("DEPT") ) 00033 { 00034 dense.lgDLWDepth = true; 00035 } 00036 else 00037 { 00038 dense.lgDLWDepth = false; 00039 } 00040 00041 p.getline(); 00042 dense.frad[0] = (realnum)p.FFmtRead(); 00043 dense.fhden[0] = (realnum)p.FFmtRead(); 00044 if( p.lgEOL() ) 00045 { 00046 fprintf( ioQQQ, " No pairs entered - can\'t interpolate.\n Sorry.\n" ); 00047 cdEXIT(EXIT_FAILURE); 00048 } 00049 00050 dense.nvals = 2; 00051 lgEnd = false; 00052 00053 /* read pairs of numbers until we find line starting with END */ 00054 /* >>chng 04 jan 27, loop to LIMTABDLAW from LIMTABD, as per 00055 * var definitions, caught by Will Henney */ 00056 while( !lgEnd && dense.nvals < LIMTABDLAW ) 00057 { 00058 p.getline(); 00059 lgEnd = p.m_lgEOF; 00060 if( !lgEnd ) 00061 { 00062 if( p.strcmp("END") == 0 ) 00063 lgEnd = true; 00064 } 00065 00066 if( !lgEnd ) 00067 { 00068 dense.frad[dense.nvals-1] = (realnum)p.FFmtRead(); 00069 dense.fhden[dense.nvals-1] = (realnum)p.FFmtRead(); 00070 dense.nvals += 1; 00071 } 00072 } 00073 --dense.nvals; 00074 00075 for( long i=1; i < dense.nvals; i++ ) 00076 { 00077 /* the radius values are assumed to be strictly increasing */ 00078 if( dense.frad[i] <= dense.frad[i-1] ) 00079 { 00080 fprintf( ioQQQ, " Radii must be in increasing order. Sorry.\n" ); 00081 cdEXIT(EXIT_FAILURE); 00082 } 00083 } 00084 } 00085 else if( p.nMatch("WIND") ) 00086 { 00087 strcpy( dense.chDenseLaw, "DLW3" ); 00088 /* This sets up a steady-state "wind" profile parametrized as in Springmann (1994): 00089 * 00090 * v(r) = v_star + (v_inf - v_0) * sqrt( Beta1 x + (1-Beta1) x^Beta2 ) 00091 * 00092 * A mass loss rate into 4pi sterradians Mdot then allows the density via continuity: 00093 * 00094 * n(r) = Mdot / ( 4Pi m_H * mu * r^2 * v(r) ) 00095 */ 00096 00097 /* The parameters must be specified in this order: 00098 * 00099 * Mdot, v_inf, Beta2, Beta1, v_0, v_star. 00100 * 00101 * Only the first three are required. The final three may be omitted right to left and 00102 * take default values Beta1 = v_0 = v_star = 0. */ 00103 00104 for( j=0; j < 6; j++ ) 00105 { 00106 dense.DensityLaw[j] = p.FFmtRead(); 00107 if( j <= 2 && p.lgEOL() ) 00108 p.NoNumb("density law element"); 00109 } 00110 } 00111 else 00112 { 00113 /* this is usual case, call dense_fabden to get density */ 00114 for( j=0; j < 10; j++ ) 00115 { 00116 dense.DensityLaw[j] = p.FFmtRead(); 00117 if( j == 0 && p.lgEOL() ) 00118 p.NoNumb("density law element"); 00119 } 00120 00121 /* set flag so we know which law to use later */ 00122 strcpy( dense.chDenseLaw, "DLW1" ); 00123 00124 /* vary option */ 00125 if( optimize.lgVarOn ) 00126 { 00127 ostringstream oss; 00128 oss << "DLAW %f" << setprecision(7); 00129 for( j=1; j < 10; j++ ) 00130 oss << " " << dense.DensityLaw[j]; 00131 strcpy( optimize.chVarFmt[optimize.nparm], oss.str().c_str() ); 00132 optimize.lgOptimizeAsLinear[optimize.nparm] = true; 00133 00134 /* index for where to write */ 00135 optimize.nvfpnt[optimize.nparm] = input.nRead; 00136 optimize.vparm[0][optimize.nparm] = (realnum)dense.DensityLaw[0]; 00137 optimize.vincr[optimize.nparm] = 0.5f; 00138 optimize.nvarxt[optimize.nparm] = 1; 00139 ++optimize.nparm; 00140 } 00141 } 00142 00143 /* set fake density to signal that density command was entered */ 00144 /* real density will be set once all input commands have been read */ 00145 /* this is necessary since density may depend on subsequent commands */ 00146 dense.SetGasPhaseDensity( ipHYDROGEN, 1.f ); 00147 00148 return; 00149 }