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 /*dense_fabden called by dlaw command, returns density for any density law */ 00004 #include "cddefines.h" 00005 #include "rfield.h" 00006 #include "dense.h" 00007 00008 /*dense_fabden called by dlaw command, returns density for any density law */ 00009 double dense_fabden(double radius, 00010 double depth) 00011 { 00012 double a , b; 00013 a = depth; 00014 b = radius; 00015 00016 fprintf(ioQQQ,"PROBLEM DISASTER The default version of dense_fabden is " 00017 "still in dense_fabden.cpp - edit and replace this file to use your version.\n"); 00018 /* return val is example - must return density */ 00019 return(a*b); 00020 00021 /* these are two examples of using this routine */ 00022 # if 0 00023 double fabden_v; 00024 double z , z0 , density , rad_midplane; 00025 DEBUG_ENTRY( "dense_fabden()" ); 00026 00027 /* this routine is called when the DLAW command is given, 00028 * and must return the hydrogen density at the current radius or depth 00029 * RADIUS is the radius, the distance from the center of symmetry. 00030 * DEPTH is the depth into the cloud, from the illuminated face 00031 * both are in cm 00032 * 00033 * radius, depth, and the array DensityLaw are double precision, although 00034 * dense_fabden itself is not 00035 * 00036 * this is one way to generate a density 00037 * dense_fabden = radius * depth 00038 * 00039 * this is how to use the parameters on the dlaw command 00040 * dense_fabden = DensityLaw(1) + radius * DensityLaw(2) 00041 * 00042 * following must be removed if this sub is to be used */ 00043 /*fprintf( ioQQQ, " the old version of dense_fabden must be deleted, and a new one put in place. Sorry.\n" );*/ 00044 00045 if( depth > radius ) 00046 TotalInsanity(); 00047 00048 rad_midplane = radius * cos( dense.DensityLaw[0] ); 00049 z = radius * sin( dense.DensityLaw[0] ); 00050 z0 = 0.047 * pow( rad_midplane/AU , 1.25 ) * AU; 00051 /* density in gm / cm3 */ 00052 density = 1.4e-9 * pow( rad_midplane/AU , -11./4. ) * sexp( POW2(z/z0) ); 00053 /* convert density to cm-3 */ 00054 fabden_v = density/(ATOMIC_MASS_UNIT * 1.2); 00055 00056 fprintf(ioQQQ, "bug dense_fabden zone %.2f radius %e density %e \n", 00057 fnzone , radius , fabden_v ); 00058 return( fabden_v ); 00059 # endif 00060 00061 #if 0 00062 00063 /* update stromgren radius to most recent value */ 00064 if( rfield.lgUSphON ) 00065 dense.DensityLaw[2] = rfield.rstrom/PARSEC; 00066 00067 /* this is the density law used in the Wen & O'Dell, 1995, ApJ 438, 784 paper */ 00068 powexp = MIN2((radius/parsec-dense.DensityLaw[2])/dense.DensityLaw[1],dense.DensityLaw[3]); 00069 fabden_v = pow(10.,dense.DensityLaw[0])*exp(powexp); 00070 00071 fabden_v = dense.DensityLaw[0]; 00072 /* just here to stop compilers from flagging unused vars */ 00073 temp = radius + depth; 00074 if( fabden_v == 0. ) 00075 { 00076 cdEXIT(EXIT_FAILURE); 00077 } 00078 else 00079 { 00080 /* when this routine is used the following branch is the correct exit */ 00081 /* following should return correct value of the density at this position, 00082 * temp is there to trick lint */ 00083 return( fabden_v*temp ); 00084 } 00085 #endif 00086 }