/*qheat1 one of Kevin Volk's quantum heating routines for grains */ #include "cddefines.h" #include "rfield.h" #include "grainvar.h" #include "phycon.h" #include "trace.h" #include "powi.h" #include "ufunct.h" #include "dintg.h" #include "nphot.h" #include "showme.h" #include "splint.h" #include "qheat1.h" void qheat1( /* this is dust type on C scale */ long int ndust1, double *tmean, double ted1) { int lgHit; long int i, i1, j, k, l; double aval[NQGRID], conv, convi, corsav[NCELL], delu[NQGRID], fract, pmax, prob[NQGRID], qtmax1, scale, sum, sum1, term, tval[NQGRID], u1[NQGRID], value, x, y; float umax, bval[NQGRID][NQGRID], umin; # ifdef DEBUG_FUN fputs( "<+>qheat1()\n", debug_fp ); # endif if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " qheat1 called grain%3ld\n", ndust1+1 ); } scale = GrainVar.qhnrat[ndust1]/GrainVar.qdep[ndust1]; for( j=0; j < NCELL; j++ ) { corsav[j] = rfield.ContBoltz[j]; } conv = 4.5911e10; /* conversion from ergs to rydberg energy units */ convi = 1./conv; /* total collisional heating in rydb/s/h atom * etot=qcu1/4. * * do j=1,nflux * value=(flux(j)+otslin(j)+outcon(j)+otscon(j))* * c dstab1(j,ndust1)*anu(j) * etot=etot+value*hden*dstAbund(j) * end do * * do j=1,nflux-1 * >>chng 97 july 17, restructure to get rid of go to */ GrainVar.qdudt0 = GrainVar.qcu1; j = 0; while( rfield.anu[j] < GrainVar.qwavec ) { /* >>chng 97 July 17, to summed continuum */ value = rfield.SummedCon[j]*GrainVar.dstab1[ndust1][j]* rfield.anu[j]; GrainVar.qdudt0 += (float)(value*phycon.hden*4.); /* end if */ ++j; } /* j now points to above qwavec */ value = rfield.SummedCon[j]*GrainVar.dstab1[ndust1][j]* rfield.anu[j]; fract = (GrainVar.qwavec - rfield.anu[j-1])/(rfield.anu[j] - rfield.anu[j-1]); value *= fract; GrainVar.qdudt0 += (float)(value*phycon.hden*4.*GrainVar.dstAbund[ndust1]); /* >>chng 97 july 17, divide out variable abundances */ /* here, qdudt0 is 1.5e-19, should have been 8.1-6 */ x = log10(convi*GrainVar.qdudt0/phycon.hden/GrainVar.dstAbund[ndust1]); splint(&GrainVar.dstems[ndust1][0],GrainVar.dsttmp,&GrainVar.dstslp[ndust1][0], NDEMS,x,&y); GrainVar.qtmin = (float)pow(10.,y); qtmax1 = MIN2(20.*ted1,GrainVar.qtmax); if( GrainVar.qtmin > GrainVar.qtmax ) { /* error condition -- set temperature to the equilibrium value and exit */ GrainVar.qtemp[0][ndust1] = (float)ted1; GrainVar.qhprob[0][ndust1] = 1.0; for( j=1; j < GrainVar.nqhg1; j++ ) { GrainVar.qtemp[j][ndust1] = (float)ted1; GrainVar.qhprob[j][ndust1] = 1.e-08f; } /* >>chng 97 mar 1, had not returned anything */ *tmean = ted1; if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " qheat1 error returns too hot\n" ); } # ifdef DEBUG_FUN fputs( " <->qheat1()\n", debug_fp ); # endif return; } x = log(qtmax1/GrainVar.qtmin)/(float)(GrainVar.nqhg1-1); i = 0; prob[0] = GrainVar.qdudt0*scale; tval[0] = GrainVar.qtmin; u1[0] = ufunct(GrainVar.qtmin,ndust1+1); for( j=1; j < GrainVar.nqhg1; j++ ) { tval[j] = GrainVar.qtmin*exp(x*(float)(j)); prob[j] = 4.*conv*dintg(tval[j],ndust1,i)*scale; u1[j] = ufunct(tval[j],ndust1+1); } delu[0] = u1[1] - u1[0]; delu[GrainVar.nqhg1-1] = u1[GrainVar.nqhg1-1] - u1[GrainVar.nqhg1-2]; bval[GrainVar.nqhg1-1][GrainVar.nqhg1-2] = (float)(prob[GrainVar.nqhg1-1] - prob[0]); for( j=1; j < (GrainVar.nqhg1 - 1); j++ ) { bval[j][j-1] = (float)(prob[j] - prob[0]); delu[j] = (u1[j+1] - u1[j-1])/2.; } for( j=0; j < (GrainVar.nqhg1 - 1); j++ ) { for( k=j+1; k < GrainVar.nqhg1; k++ ) { umin = (float)(u1[k] - delu[k]/2. - u1[j]); umax = (float)(u1[k] + delu[k]/2. - u1[j]); if( (k+1) == GrainVar.nqhg1 ) umax = 100000.; nphot(umin,umax,&bval[j][k],ndust1); } } for( k=0; k < GrainVar.nqhg1; k++ ) { term = log10(tval[k]); /* interpolate to find the energy loss value at t = tval(k) * */ i1 = 0; lgHit = FALSE; j = 1; while( !lgHit && j < NDEMS ) { /* >>chng 97 jan 3, by gjf, to get rid of statement number * do j=1,NDEMS-1 */ if( term > GrainVar.dsttmp[j-1] && term <= GrainVar.dsttmp[j] ) { fract = (term - GrainVar.dsttmp[j-1])/(GrainVar.dsttmp[j] - GrainVar.dsttmp[j-1]); value = GrainVar.dstems[ndust1][j-1] + fract*(GrainVar.dstems[ndust1][j] - GrainVar.dstems[ndust1][j-1]); i1 = 1; lgHit = TRUE; } j += 1; } if( i1 == 1 ) { aval[k] = scale*conv*(pow(10.,value)); } else { aval[k] = scale*conv*(pow(10.,GrainVar.dstems[ndust1][0])); } } for( j=GrainVar.nqhg1-1; j >= 1; j-- ) { aval[j] -= aval[0]; aval[j] = delu[j]/aval[j]; } aval[0] = 0.; for( j=0; j < (GrainVar.nqhg1 - 2); j++ ) { for( k=GrainVar.nqhg1 - 2; k > j; k-- ) /* >>>99nov14 logic error in below? */ /*for( k=GrainVar.nqhg1 - 2; k >= (j); k-- )*/ { bval[j][k] += bval[j][k+1]; } } prob[0] = 1.0; pmax = 1.0; sum = prob[0]; for( k=1; k < GrainVar.nqhg1; k++ ) { prob[k] = 0.; for( j=0; j < k; j++ ) { prob[k] += prob[j]*bval[j][k]; } prob[k] *= aval[k]; sum += prob[k]; if( prob[k] > pmax ) pmax = prob[k]; if( pmax > 1.e06 ) { value = 1./pmax; for( l=0; l < k+1; l++ ) { prob[l] *= value; } sum *= value; pmax = 1.0; } } sum = 1./sum; for( j=0; j < GrainVar.nqhg1; j++ ) { prob[j] *= sum; } sum = 0.; sum1 = 0.; for( j=0; j < GrainVar.nqhg1; j++ ) { sum += prob[j]; sum1 += prob[j]*(powi(tval[j],4)); } /* this is the final temperature */ *tmean = pow(sum1,0.25); /* check that insanity has not happened */ if( *tmean <= 0 ) { fprintf(ioQQQ, " qheat1 insanity, non-positive grain temperature. This is impossible.\n"); ShowMe(); puts( "[Stop in qheat1]" ); exit(1); } for( j=0; j < GrainVar.nqhg1; j++ ) { GrainVar.qtemp[j][ndust1] = (float)tval[j]; GrainVar.qhprob[j][ndust1] = (float)prob[j]; /* debugging aid * write(42,398) qtemp(ndust1,j),qhprob(ndust1,j) * 398 format(1x,0pf8.3,1x,1pe12.5) */ } for( j=0; j < NCELL; j++ ) { rfield.ContBoltz[j] = (float)corsav[j]; } if( trace.lgTrace && trace.lgDustBug ) { fprintf( ioQQQ, " qheat1 returns%10.2e%10.2e\n", *tmean, ted1 ); } # ifdef DEBUG_FUN fputs( " <->qheat1()\n", debug_fp ); # endif return; }