/*dcharg compute grains charge */ #include #include "cddefines.h" #include "sign.h" #define DMAXX 0.05 #define TOLER 0.005 #include "physconst.h" #include "grainvar.h" #include "trace.h" #include "converge.h" #include "ipoint.h" #include "grngam.h" #include "dcharg.h" /*drecom compute electron recombination onto grains surface */ double drecom( /* nd is pointer to the grain species */ long int nd , /* this is collider rate for neutral species */ double *RecNeutral); void dcharg(void) { int lgBigError; long int ip , lopbig, luplim, nd; double change, dnew, error, phtion, recom, RecNeutral , rerror, thres, xxx; static double oldpht[NDUST], oldpot[NDUST], oldrec[NDUST], slope[NDUST]; # ifdef DEBUG_FUN fputs( "<+>dcharg()\n", debug_fp ); # endif if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " dcharg called lgSearch%2c\n", TorF(conv.lgSearch) ); } for( nd=0; nd < NDUST; nd++ ) { if( GrainVar.lgDustOn1[nd] ) { if( conv.lgSearch ) { if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " dchar search %10.10s\n", chGrainVar.chDstLab[nd] ); } luplim = 100; GrainVar.dstpot[nd] = 0.; /* DustWorkFcn is neutral surface work function in ryd * DSTPOT is floating grain pot in Ryd */ thres = MAX2(GrainVar.DustWorkFcn[nd],GrainVar.DustWorkFcn[nd]+ GrainVar.dstpot[nd]); /* ipDwork is the actualy pointer to threshold generated * by calling ipContSafe in setup - it could be a few cells * to either side to the correct cell if several grains * share the same ip, making max necessary below */ ip = ipoint(thres); GrainVar.ipDustPot[nd] = MAX2(GrainVar.ipDwork[nd] , ip ); oldpot[nd] = GrainVar.dstpot[nd]; /* functiong rngam returns photoionization rate for this grain species*/ oldpht[nd] = grngam(nd); /* the function drecom returns the recombination rate * for this grain and charge */ oldrec[nd] = drecom(nd, &RecNeutral ); GrainVar.dstpot[nd] += 0.05f; } else { /* this is limit to number of loops we will try before giving up on * getting the grain charge right, set to 100 when in search phase */ luplim = 10; } /* find dust charge; do sum over photon number */ lopbig = 0; lgBigError = TRUE; phtion = -DBL_MAX; recom = -DBL_MAX; while( lopbig <= luplim && lgBigError ) { lopbig += 1; thres = MAX2(GrainVar.DustWorkFcn[nd],GrainVar.DustWorkFcn[nd]+ GrainVar.dstpot[nd]); GrainVar.ipDustPot[nd] = ipoint(thres); /* functiong rngam returns photoionization rate for this grain species*/ phtion = grngam(nd); /* the function drecom returns the recombination rate * for this grain and charge */ recom = drecom(nd, &RecNeutral ); /* this is the difference between photoionization and recombination * rates - we want this to be small */ error = phtion - recom; /* >>chng 97 feb 28, crashed when both phtion and recom zero - which *happens when no ionizing radiation */ /* >>chng 99 jun 16, grainlte.in has no ionizing rad at all, and loop * never converged. add second test for convergence, since first test * always failts when phtion is zero and recom is relative to recom */ if( phtion > SMALLFLOAT ) { rerror = fabs(error) / (0.5*MAX2(SMALLFLOAT,phtion+recom)); } else { /* RecNeutral evaluated in drecom (which is below) and is collision * rate were the grains neutral */ rerror = fabs(error) / MAX2(SMALLFLOAT , RecNeutral); } if( rerror > TOLER ) { if( GrainVar.dstpot[nd] != oldpot[nd] ) { /* only re-estimate slope if new continuum * >>chng 96 apr 26, as per Peter van Hoof comment that div by zero can occur */ xxx = (phtion - oldpht[nd]) - (recom - oldrec[nd]); if( xxx != 0. ) slope[nd] = xxx/(GrainVar.dstpot[nd] - oldpot[nd]); } oldpot[nd] = GrainVar.dstpot[nd]; oldpht[nd] = phtion; oldrec[nd] = recom; /* here is reason for above xxx, do not want zero slope: * really the change should be here not above?? */ dnew = -error/slope[nd]; change = fabs(dnew); change = MIN2(change,DMAXX*thres); change = sign(change,dnew); GrainVar.dstpot[nd] += (float)change; lgBigError = TRUE; } else { oldpot[nd] = GrainVar.dstpot[nd]; oldpht[nd] = phtion; oldrec[nd] = recom; lgBigError = FALSE; } if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " dcharg loop %3ld%10.10s pot%10.2e error%10.2e Pho rec Rate%9.1e%9.1e\n", lopbig, chGrainVar.chDstLab[nd], GrainVar.dstpot[nd], rerror, phtion, recom ); } } if( lgBigError && !conv.lgSearch ) { fprintf( ioQQQ," dcharg did not converge ionization of grain species %li\n", nd); /* increment the counter */ ++conv.nGrainFail; } GrainVar.dstq[nd] = (float)((GrainVar.dstpot[nd]*EVRYD)/0.0048); if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " Grain potential:%10.10s%10.2eRyd, gamma:%10.2e recom:%10.2e Threshold:%10.2eRyd\n", chGrainVar.chDstLab[nd], GrainVar.dstpot[nd], phtion, recom, GrainVar.DustWorkFcn[nd] + GrainVar.dstpot[nd] ); } } } if( trace.lgTrace && (!trace.lgDustBug) ) { fprintf( ioQQQ, " Grain potentials:" ); for( nd=0; nd < NDUST; nd++ ) { fprintf( ioQQQ, "%10.2e", GrainVar.dstpot[nd] ); } fprintf( ioQQQ, "\n" ); } # ifdef DEBUG_FUN fputs( " <->dcharg()\n", debug_fp ); # endif return; } /*drecom compute electron recombination onto grains surface */ #define STICK 1. #include "phycon.h" #include "sexp.h" #include "tepowers.h" double drecom( /* nd is pointer to the grain species */ long int nd, /* this is collider rate for neutral species */ double *RecNeutral) { double drecom_v, eta, psi, ve; # ifdef DEBUG_FUN fputs( "<+>drecom()\n", debug_fp ); # endif /* electron sticking prob */ /* coef is qe(esu) / 300 / k to convert potential in volts to PSI * in eV here * this is defined for impact by electrons */ psi = -1.*(GrainVar.dstpot[nd]*EVRYD)*1.1606e4/phycon.te; if( psi <= 0. ) { eta = 1. - psi; } else { eta = sexp(psi); } /* this is rec coef for neutral grain; corr fac next * VE is mean electron velocity */ ve = TePowers.sqrte*6.2124e5; /* this is neutral surface recomb rate for electrons */ *RecNeutral = ve*phycon.eden*STICK; /* now take into account screening by grain charge */ drecom_v = *RecNeutral*eta; # ifdef DEBUG_FUN fputs( " <->drecom()\n", debug_fp ); # endif return( drecom_v ); }