cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ParseRadius parse the radius command */ 00004 #include "cddefines.h" 00005 /*#define PARSCL 18.489396*/ 00006 #include "physconst.h" 00007 #include "optimize.h" 00008 #include "radius.h" 00009 #include "iterations.h" 00010 #include "input.h" 00011 #include "parse.h" 00012 00013 void ParseRadius(char *chCard, 00014 realnum *ar1) 00015 { 00016 bool lgEOL, 00017 lgR2set, 00018 lgRLog; 00019 long int i; 00020 double a, 00021 convl; 00022 00023 DEBUG_ENTRY( "ParseRadius()" ); 00024 00025 /* log of inner and outer radii, default second=infinity, 00026 * if R2<R1 then R2=R1+R2 00027 * there is an optional keyword, "PARSEC" on the line, to use PC as units */ 00028 if( nMatch("PARS",chCard) ) 00029 { 00030 /*>>chng 06 mar 18, from above to below */ 00031 /*convl = PARSCL;*/ 00032 convl = log10( PARSEC ); 00033 } 00034 else 00035 { 00036 convl = 0.; 00037 } 00038 00039 /* if linear appears on line, then radius is linear, otherwise, log */ 00040 if( nMatch("LINE",chCard) ) 00041 { 00042 lgRLog = false; 00043 } 00044 else 00045 { 00046 lgRLog = true; 00047 } 00048 00049 i = 5; 00050 *ar1 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00051 if( lgEOL ) 00052 { 00053 fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" ); 00054 cdEXIT(EXIT_FAILURE); 00055 } 00056 00057 /* option for linear or log radius */ 00058 if( lgRLog ) 00059 { 00060 *ar1 += (realnum)convl; 00061 } 00062 else 00063 { 00064 if( *ar1 > 0. ) 00065 { 00066 *ar1 = (realnum)(log10(*ar1) + convl); 00067 } 00068 else 00069 { 00070 fprintf(ioQQQ,"The first radius is negative and linear radius is set - this is impossible.\n"); 00071 cdEXIT(EXIT_FAILURE); 00072 } 00073 } 00074 00075 if( *ar1 > 37 || *ar1 < -37. ) 00076 { 00077 fprintf(ioQQQ,"WARNING - the log of the radius is %e - this is too big and I would crash.\n", *ar1 ); 00078 fprintf(ioQQQ," Sorry.\n" ); 00079 cdEXIT(EXIT_FAILURE); 00080 } 00081 00082 radius.Radius = pow((realnum)10.f ,*ar1); 00083 radius.lgRadiusKnown = true; 00084 00085 /* check for second number, which indicates thickness or outer radius of model */ 00086 a = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00087 if( lgEOL ) 00088 { 00089 /* not set */ 00090 lgR2set = false; 00091 } 00092 else 00093 { 00094 /* outer radius is set, */ 00095 lgR2set = true; 00096 00097 /* log or linear option is still in place */ 00098 if( lgRLog ) 00099 { 00100 a += convl; 00101 } 00102 00103 else 00104 { 00105 /* linear radius - convert to log but first make sure that a is > 0 */ 00106 if( a > 0. ) 00107 { 00108 a = log10(a) + convl; 00109 } 00110 else 00111 { 00112 fprintf(ioQQQ,"The second radius is negative and linear radius is set - this is impossible.\n"); 00113 cdEXIT(EXIT_FAILURE); 00114 } 00115 } 00116 00117 if( a > 37 || a < -37. ) 00118 { 00119 fprintf(ioQQQ,"WARNING - the log of the outer radius is %e - this is too big and I will soon crash.\n", a ); 00120 /* flush buffers since we shall soon throw an fpe */ 00121 fflush( ioQQQ ); 00122 } 00123 a = pow(10.,a); 00124 /* check whether it was thickness or outer radius, 00125 * we want router to be total thickness of modeled region, 00126 * NOT outer radius */ 00127 if( a > radius.Radius ) 00128 { 00129 radius.router[0] = a - radius.Radius; 00130 } 00131 00132 else 00133 { 00134 radius.router[0] = a; 00135 } 00136 00137 for( i=1; i < iterations.iter_malloc; i++ ) 00138 { 00139 radius.router[i] = radius.router[0]; 00140 } 00141 } 00142 00143 /* vary option */ 00144 if( optimize.lgVarOn ) 00145 { 00146 /* pointer to where to write */ 00147 optimize.nvfpnt[optimize.nparm] = input.nRead; 00148 optimize.vincr[optimize.nparm] = 0.5; 00149 00150 /* flag saying second outer radius or thickness was set */ 00151 if( lgR2set ) 00152 { 00153 strcpy( optimize.chVarFmt[optimize.nparm], "RADIUS %f depth or outer R %f" ); 00154 optimize.nvarxt[optimize.nparm] = 2; 00155 /* second number is thickness or outer radius */ 00156 optimize.vparm[1][optimize.nparm] = (realnum)log10(radius.router[0]); 00157 fprintf(ioQQQ, 00158 " WARNING - outer radius or thickness was set with a variable radius.\n"); 00159 fprintf(ioQQQ, 00160 " The interpretation of the second number can change from radius to depth as radius changes.\n"); 00161 fprintf(ioQQQ, 00162 " Do not use the second parameter unless you are quite certain that you know what you are doing.\n"); 00163 fprintf(ioQQQ, 00164 " Consider using the STOP THICKNESS command instead.\n"); 00165 } 00166 else 00167 { 00168 strcpy( optimize.chVarFmt[optimize.nparm], "RADIUS= %f" ); 00169 optimize.nvarxt[optimize.nparm] = 1; 00170 } 00171 00172 /* log of radius is first number */ 00173 optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.Radius); 00174 ++optimize.nparm; 00175 } 00176 return; 00177 }