/*SumContinuum sum flux, otscon, otslin, outcon, outlin, to form SummeDif, SummedCon SummedOcc */ #include "cddefines.h" #include "rfield.h" #include "opac.h" #include "sumcontinuum.h" void SumContinuum(double *SumOTS , double *RatioOTSInc , long int * ipOTSchange) { long int i; /* when setting new to sum of current and old, this is the fraction of old to take*/ float FracOld , FracNew , changeOTS , OpacRatio ; # ifdef DEBUG_FUN fputs( "<+>SumContinuum()\n", debug_fp ); # endif *SumOTS = 0.; *RatioOTSInc = 0.; # define FRACOLD 0.25f /* fraction old and new ots rates */ FracOld = FRACOLD; FracNew = 1.f - FracOld ; /* remember largest change in ots rates */ changeOTS = 0.; *ipOTSchange = -100000; /*fprintf(ioQQQ," sumcontinuum enter 589=%e\n", rfield.otslin[589]);*/ for( i=0; i < rfield.nflux; i++ ) { /* remember which energy had largest change in ots rates */ if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS) { *ipOTSchange = i; changeOTS = (float) fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i]); } /* ots rates were evaluated with previous opacities, now have different * opacities since AddOpac just called. Rescale OTS rates to current opacity. * tests with pure h, he constant temperature clouds shows this decreases*/ /* >>chng 99 apr 9, rescale ots */ OpacRatio = 1.; OpacRatio = (float)(opac.OldOpacSave[i]/opac.opac[i]); /* the New vectors are the ones updated by the AddOTS routines, * and have the current OTS rates. The otscon and otslin vectors * have the rates from the previous refresh of the New vectors*/ /* here is the refresh. all were initially set to zero in call to * ZeroOTSContinuum (below) from zero */ rfield.otscon[i] = (rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew)*OpacRatio; /* zero out the new otscon vector*/ rfield.otsconNew[i] = 0.; /* current ots will be average of old and new */ rfield.otslin[i] = (rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew)*OpacRatio; /* zero out the new otslin vector*/ rfield.otslinNew[i] = 0.; /* remember sum of ots rates for convergence criteria */ *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opac[i]; /* also sum incident continuum for comparison */ *RatioOTSInc += rfield.flux[i]*opac.opac[i]; rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outcon[i] + rfield.outlin[i]; rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i]; } /* form ratio of ots interaction rates to incident continuum interaction, * shows relative importance of ots */ if( *RatioOTSInc > 0. ) { *RatioOTSInc = *SumOTS / *RatioOTSInc ; } else { *RatioOTSInc = 0.; } /*fprintf(ioQQQ," sumcontinuum exit 589=%e\n", rfield.otslin[589]);*/ /*lint -e506 constant boolean */ /*lint -e774 boolenan with constant if */ { /* following should be set true to print strongest ots contributors */ enum {DEBUG = FALSE}; if( DEBUG /*&& (IterCnt.iter==2)*/ ) { double BigOTS; int ipOTS=0; BigOTS = -1.; for( i=0; i < rfield.nflux; i++ ) { if( (rfield.otscon[i] + rfield.otslin[i])*opac.opac[i] > BigOTS ) { BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opac[i]; ipOTS = i; } } fprintf(ioQQQ," sumcontinuunots %10.2e %10.2e %10.2e\n", rfield.anu[ipOTS] , opac.opac[ipOTS],BigOTS ); } } /*lint +e506 constant boolean */ /*lint +e774 boolenan with constant if */ # ifdef DEBUG_FUN fputs( " <->SumContinuum()\n", debug_fp ); # endif return; } /*ZeroOTSContinuum - zero out some vectors - * this is only called when code initialized by setcon */ void ZeroOTSContinuum( void ) { long int i; # ifdef DEBUG_FUN fputs( "<+>ZeroOTSContinuum()\n", debug_fp ); # endif /* this loop goes up to nflux itself (<=) since the highest cell * will be used to pass unity through the code to verify integrations */ for( i=0; i <= rfield.nflux; i++ ) { rfield.SummedDif[i] = 0.; /* the four main ots vectors */ rfield.otscon[i] = 0.; rfield.otslin[i] = 0.; rfield.otslinNew[i] = 0.; rfield.otsconNew[i] = 0.; rfield.outcon[i] = 0.; rfield.outlin[i] = 0.; rfield.SummedDif[i] = 0.; /* "zero" for the summed con will be just the incident radiation field */ rfield.SummedCon[i] = rfield.flux[i]; rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i]; } # ifdef DEBUG_FUN fputs( " <->ZeroOTSContinuum()\n", debug_fp ); # endif return; }