/*tabden interpolate on table of points for density with dlaw table command, by K Volk */ #include #include "cddefines.h" #include "abuntablfac.h" #include "tabden.h" double tabden(double r0, double depth) { int lgHit; long int j; double frac, tabden_v, x; # ifdef DEBUG_FUN fputs( "<+>tabden()\n", debug_fp ); # endif /*interpolate on table of points for density with dlaw table command, by K Volk *each line is log radius and H density per cc. */ /*begin sanity check */ if( r0 <= 0. || depth <= 0. ) { fprintf( ioQQQ, " tabden called with insane depth, radius, =%10.2e%10.2e\n", depth, r0 ); } /*end sanity check */ /* interpolate on radius or depth? */ if( AbunTablFac.lgDLWDepth ) { /* depth key appeared = we want depth */ x = log10(depth); } else { /* use radius */ x = log10(r0); } /* set to impossible value, will crash if not reset */ tabden_v = -DBL_MAX; if( x < AbunTablFac.frad[0] || x >= AbunTablFac.frad[AbunTablFac.nvals-1] ) { fprintf( ioQQQ, " requested radius outside range of tabden\n" ); fprintf( ioQQQ, " radius was%10.2e min, max=%10.2e%10.2e\n", x, AbunTablFac.frad[0], AbunTablFac.frad[AbunTablFac.nvals-1] ); puts( "[Stop in tabden]" ); exit(1); } else { lgHit = FALSE; j = 1; while( !lgHit && j <= AbunTablFac.nvals - 1 ) { if( AbunTablFac.frad[j-1] <= x && AbunTablFac.frad[j] > x ) { frac = (x - AbunTablFac.frad[j-1])/(AbunTablFac.frad[j] - AbunTablFac.frad[j-1]); tabden_v = AbunTablFac.fhden[j-1] + frac*(AbunTablFac.fhden[j] - AbunTablFac.fhden[j-1]); lgHit = TRUE; } j += 1; } if( !lgHit ) { fprintf( ioQQQ, " radius outran dlaw table scale, requested=%6.2f largest=%6.2f\n", x, AbunTablFac.frad[AbunTablFac.nvals-1] ); puts( "[Stop in tabden]" ); exit(1); } } /* got it, now return value, not log of density */ tabden_v = pow(10.,tabden_v); # ifdef DEBUG_FUN fputs( " <->tabden()\n", debug_fp ); # endif return( tabden_v ); }