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_tabden interpolate on table of points for density with dlaw table command, by K Volk */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 00007 double dense_tabden(double r0, 00008 double depth) 00009 { 00010 bool lgHit; 00011 long int j; 00012 double frac, 00013 tabden_v, 00014 x; 00015 00016 DEBUG_ENTRY( "dense_tabden()" ); 00017 /*interpolate on table of points for density with dlaw table command, by K Volk 00018 *each line is log radius and H density per cc. */ 00019 00020 /*begin sanity check */ 00021 if( r0 <= 0. || depth <= 0. ) 00022 { 00023 fprintf( ioQQQ, " dense_tabden called with insane depth, radius, =%10.2e%10.2e\n", 00024 depth, r0 ); 00025 } 00026 /*end sanity check */ 00027 00028 /* interpolate on radius or depth? */ 00029 if( dense.lgDLWDepth ) 00030 { 00031 /* depth key appeared = we want depth */ 00032 x = log10(depth); 00033 } 00034 else 00035 { 00036 /* use radius */ 00037 x = log10(r0); 00038 } 00039 00040 /* set to impossible value, will crash if not reset */ 00041 tabden_v = -DBL_MAX; 00042 00043 if( x < dense.frad[0] || x >= dense.frad[dense.nvals-1] ) 00044 { 00045 fprintf( ioQQQ, " requested radius outside range of dense_tabden\n" ); 00046 fprintf( ioQQQ, " radius was%10.2e min, max=%10.2e%10.2e\n", 00047 x, dense.frad[0], dense.frad[dense.nvals-1] ); 00048 cdEXIT(EXIT_FAILURE); 00049 } 00050 else 00051 { 00052 lgHit = false; 00053 j = 1; 00054 00055 while( !lgHit && j <= dense.nvals - 1 ) 00056 { 00057 if( dense.frad[j-1] <= (realnum)x && dense.frad[j] > (realnum)x ) 00058 { 00059 frac = (x - dense.frad[j-1])/(dense.frad[j] - 00060 dense.frad[j-1]); 00061 tabden_v = dense.fhden[j-1] + frac*(dense.fhden[j] - 00062 dense.fhden[j-1]); 00063 lgHit = true; 00064 } 00065 j += 1; 00066 } 00067 00068 if( !lgHit ) 00069 { 00070 fprintf( ioQQQ, " radius outran dlaw table scale, requested=%6.2f largest=%6.2f\n", 00071 x, dense.frad[dense.nvals-1] ); 00072 cdEXIT(EXIT_FAILURE); 00073 } 00074 } 00075 00076 /* got it, now return value, not log of density */ 00077 tabden_v = pow(10.,tabden_v); 00078 return( tabden_v ); 00079 }