00001
00002
00003
00004 #include "cddefines.h"
00005 #include "input.h"
00006 #include "dense.h"
00007 #include "optimize.h"
00008 #include "parser.h"
00009 #include "cosmology.h"
00010
00011 void ParseHDEN(Parser &p )
00012 {
00013 DEBUG_ENTRY( "ParseHDEN()" );
00014
00015 if( dense.gas_phase[ipHYDROGEN] > 0. )
00016 {
00017 fprintf( ioQQQ, " PROBLEM DISASTER More than one density command was entered.\n" );
00018 cdEXIT(EXIT_FAILURE);
00019 }
00020
00021
00022 dense.gas_phase[ipHYDROGEN] = (realnum)p.FFmtRead();
00023 if( p.lgEOL() )
00024 {
00025 fprintf( ioQQQ, " DISASTER The density MUST be entered with this command. STOP\n" );
00026 cdEXIT(EXIT_FAILURE);
00027 }
00028
00029
00030 if( ! p.nMatch("LINE") )
00031 {
00032
00033 if( dense.gas_phase[ipHYDROGEN] > log10(FLT_MAX) ||
00034 dense.gas_phase[ipHYDROGEN] < log10(FLT_MIN) )
00035 {
00036 fprintf(ioQQQ,
00037 " DISASTER - the log of the entered hydrogen density is %.3f - too extreme for this processor.\n",
00038 dense.gas_phase[ipHYDROGEN]);
00039 if( dense.gas_phase[ipHYDROGEN] > 0. )
00040 fprintf(ioQQQ,
00041 " DISASTER - the log of the largest hydrogen density this processor can do is %.3f.\n",
00042 log10(FLT_MAX) );
00043 else
00044 fprintf(ioQQQ,
00045 " DISASTER - the log of the smallest hydrogen density this processor can do is %.3f.\n",
00046 log10(FLT_MIN) );
00047 fprintf(ioQQQ," Sorry.\n" );
00048
00049 cdEXIT(EXIT_FAILURE);
00050 }
00051
00052 dense.gas_phase[ipHYDROGEN] = (realnum)pow((realnum)10.f,dense.gas_phase[ipHYDROGEN]);
00053 }
00054
00055 if( dense.gas_phase[ipHYDROGEN] > MAX_DENSITY )
00056 {
00057 fprintf( ioQQQ, "This density is too high. This version of Cloudy does not permit densities greater than %e cm-3.\n", MAX_DENSITY );
00058 cdEXIT(EXIT_FAILURE);
00059 }
00060
00061 if( dense.gas_phase[ipHYDROGEN] <= 0. )
00062 {
00063 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density must be > 0.\n" );
00064 cdEXIT(EXIT_FAILURE);
00065 }
00066
00067
00068 dense.den0 = dense.gas_phase[ipHYDROGEN];
00069
00070
00071 dense.DensityPower = (realnum)p.FFmtRead();
00072
00073 if( !p.lgEOL() )
00074 {
00075
00076
00077 if( p.nMatch("COLU") )
00078 {
00079
00080
00081
00082 dense.rscale = (realnum)pow(10.,p.FFmtRead());
00083 if( p.lgEOL() )
00084 {
00085 fprintf( ioQQQ, " The column density MUST be set if the col den option is to be used.\n" );
00086 cdEXIT(EXIT_FAILURE);
00087 }
00088 strcpy( dense.chDenseLaw, "POWC" );
00089 }
00090 else if( p.nMatch("DEPT") )
00091 {
00092
00093 dense.rscale = (realnum)pow(10.,p.FFmtRead());
00094 if( p.lgEOL() )
00095 {
00096 fprintf( ioQQQ, " The scale depth MUST be set if the depth option is to be used.\n" );
00097 cdEXIT(EXIT_FAILURE);
00098 }
00099 strcpy( dense.chDenseLaw, "POWD" );
00100 }
00101 else
00102 {
00103
00104 strcpy( dense.chDenseLaw, "POWR" );
00105 }
00106 }
00107
00108
00109 if( optimize.lgVarOn )
00110 {
00111
00112 optimize.nvfpnt[optimize.nparm] = input.nRead;
00113 optimize.vparm[0][optimize.nparm] = (realnum)log10(dense.gas_phase[ipHYDROGEN]);
00114 optimize.vincr[optimize.nparm] = 1.;
00115
00116
00117
00118 if( strcmp(dense.chDenseLaw ,"CDEN") == 0 ||
00119 strcmp(dense.chDenseLaw ,"CPRE") == 0 ||
00120 strcmp(dense.chDenseLaw ,"WIND") == 0 )
00121 {
00122 strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f LOG" );
00123 optimize.nvarxt[optimize.nparm] = 1;
00124 }
00125
00126
00127 else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
00128 {
00129 strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f LOG, power=%f" );
00130 optimize.nvarxt[optimize.nparm] = 2;
00131 optimize.vparm[1][optimize.nparm] = dense.DensityPower;
00132 }
00133
00134
00135 else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
00136 {
00137 strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f LOG, power=%f, column=%f" );
00138 optimize.nvarxt[optimize.nparm] = 3;
00139 optimize.vparm[1][optimize.nparm] = dense.DensityPower;
00140 optimize.vparm[2][optimize.nparm] = (realnum)log10(dense.rscale);
00141 }
00142
00143
00144 else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
00145 {
00146 strcpy( optimize.chVarFmt[optimize.nparm], "HDEN=%f LOG, power=%f, depth=%f" );
00147 optimize.nvarxt[optimize.nparm] = 3;
00148 optimize.vparm[1][optimize.nparm] = dense.DensityPower;
00149 optimize.vparm[2][optimize.nparm] = (realnum)log10(dense.rscale);
00150 }
00151
00152
00153 else
00154 {
00155 fprintf( ioQQQ, " Internal error in HDEN\n" );
00156 cdEXIT(EXIT_FAILURE);
00157 }
00158 ++optimize.nparm;
00159 }
00160 return;
00161 }