46 static bool lgPntEval =
false;
57 fprintf(
ioQQQ,
" ContCreateMesh called, not evaluating.\n" );
60 for( i=0; i <
rfield.nupper; i++ )
72 fprintf(
ioQQQ,
" ContCreateMesh called first time.\n" );
125 fprintf(
ioQQQ,
" Currently the arrays that hold interpolated tables can only hold %i points.\n",
NCELL);
126 fprintf(
ioQQQ,
" This continuum mesh really needs to have %li points.\n",
rfield.nupper);
127 fprintf(
ioQQQ,
" Please increase the value of NCELL in rfield.h and recompile.\n Sorry.");
150 for( i=0;i<(
geometry.nend_max+1); ++i )
155 for( i=0;i<(
geometry.nend_max+1); ++i )
157 for(
long j=0; j<
rfield.nupper; ++j)
159 rfield.ConSourceFcnLocal[i][j] = 1.;
200 for( i=1; i<
rfield.nupper-1; ++i )
208 for( i=0; i <
rfield.nupper; i++ )
216 alf = 1./(1. +
rfield.anu[i]*(1.1792e-4 + 7.084e-10*
rfield.anu[i]));
217 bet = 1. - alf*
rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*
rfield.anu[i])/4.;
227 rfield.ipnt_coarse_2_fine[i] = 0;
234 rfield.ipnt_coarse_2_fine[i] = 0;
244 rfield.ipnt_coarse_2_fine[i] = ipnt;
250 rfield.resetCoarseTransCoef();
276 ASSERT( fenlo>0. && fenhi>0. && resolv>0. );
279 nbin = (
long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);
288 if( *ipnt > 0 && fabs(1.-fenlo/
continuum.filbnd[*ipnt]) > 1e-4 )
290 fprintf(
ioQQQ,
" FILL improper bounds.\n" );
291 fprintf(
ioQQQ,
" ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n",
316 if( (*n0 + nbin-2) >
rfield.nupper )
318 fprintf(
ioQQQ,
" Fill would need %ld cells to get to an energy of %.3e\n",
320 fprintf(
ioQQQ,
" This is a major logical error in fill.\n");
326 for( i=0; i < nbin; i++ )
329 aaa = pow( 10. , bbb );
343 " FILL range%2ld from%10.3e to%10.3eR in%4ld cell; ener res=%10.3e WIDTOT=%10.3e\n",
351 fprintf(
ioQQQ,
" The requested range was%10.3e%10.3e The requested resolution was%10.3e\n",
352 fenlo, fenhi, resolv );
456 rfield.fine_ener_hi = 1500.f;
474 double TeLowestFineOpacity = 1e4;
475 rfield.fine_opac_velocity_width =
483 rfield.ipFineConVelShift = 0;
489 rfield.nfine_malloc = (long)(log10(
rfield.fine_ener_hi /
rfield.fine_ener_lo ) / log10( 1. +
rfield.fine_resol ) );
490 if(
rfield.nfine_malloc <= 0 )
506 for( i=0;i<
rfield.nfine_malloc; ++i )
514 for( i=0; i<
rfield.nupper; ++i)
552 rfield.ipnt_coarse_2_fine = (
long int*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
long int)) );
585 for( i = 1; i <=
LIMELM; i++ )
599 for( i=0; i<
rfield.nupper; ++i)
610 for( i=0; i<
rfield.nupper; ++i)
627 opac.opacity_abs = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
628 opac.opacity_abs_savzon1 = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
629 opac.OldOpacSave = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
630 opac.opacity_sct = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
632 opac.opacity_sct_savzon1 = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
634 opac.FreeFreeOpacity = (
double*)
MALLOC((
size_t)(
rfield.nupper*
sizeof(
double)) );
669 fprintf(
ioQQQ,
" read_continuum_mesh opening continuum_mesh.ini:");
671 ioDATA =
open_data(
"continuum_mesh.ini",
"r" );
676 fprintf(
ioQQQ,
" read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
682 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
686 if( chLine[0] !=
'#')
701 ((
double *)
MALLOC( (
size_t)(
continuum.nStoredBands+1)*
sizeof(
double )));
703 ((
double *)
MALLOC( (
size_t)(
continuum.nStoredBands+1)*
sizeof(
double )));
706 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
708 fprintf(
ioQQQ,
" read_continuum_mesh could not rewind continuum_mesh.ini.\n");
715 fprintf(
ioQQQ,
" read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
721 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
722 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
723 i3 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
728 if( i1 == 1 && i2 == 9 && i3 == 29 )
732 else if( i1 == 10 && i2 == 8 && i3 == 8 )
739 " read_continuum_mesh: the version of continuum_mesh.ini is not supported.\n" );
741 " I found version number %li %li %li.\n" ,
743 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
750 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
753 if( chLine[0] !=
'#')
764 fprintf(
ioQQQ,
"DISASTER PROBLEM continuum_mesh.ini has a non-positive number.\n");
783 for( i=1; i<
continuum.nStoredBands-1; ++i )
788 " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
795 " read_continuum_mesh: The last continuum array energies must be zero.\n" );
815 unsigned long n=(
unsigned long)(ihi-lo+1);
816 memset(&
rfield.OccNumbDiffCont[lo] , 0 , n*
sizeof(
realnum) );
817 memset(&
rfield.OccNumbContEmitOut[lo] , 0 , n*
sizeof(
realnum) );
818 memset(&
rfield.ContBoltz[lo] , 0 , n*
sizeof(
double) );
822 memset(&
rfield.ConEmitReflec[0][lo] , 0 , n*
sizeof(
realnum) );
826 memset(&
rfield.SummedCon[lo] , 0 , n*
sizeof(
double) );
827 memset(&
rfield.OccNumbBremsCont[lo] , 0 , n*
sizeof(
realnum) );
830 memset(&
rfield.flux_total_incident[0][lo] , 0 , n*
sizeof(
realnum) );
831 memset(&
rfield.flux_beam_const_save[lo] , 0 , n*
sizeof(
realnum) );
832 memset(&
rfield.flux_time_beam_save[lo] , 0 , n*
sizeof(
realnum) );
833 memset(&
rfield.flux_isotropic_save[lo] , 0 , n*
sizeof(
realnum) );
842 memset(&
rfield.ConOTS_local_OTS_rate[lo], 0 , n*
sizeof(
realnum) );
843 memset(&
rfield.ConOTS_local_photons[lo] , 0 , n*
sizeof(
realnum) );
844 memset(&
opac.OldOpacSave[lo] , 0 , n*
sizeof(
double) );
845 memset(&
opac.opacity_abs[lo] , 0 , n*
sizeof(
double) );
846 memset(&
opac.opacity_sct[lo] , 0 , n*
sizeof(
double) );
847 memset(&
opac.albedo[lo] , 0 , n*
sizeof(
double) );
848 memset(&
opac.FreeFreeOpacity[lo] , 0 , n*
sizeof(
double) );
851 memset( &
opac.E2TauAbsTotal[lo] , 0 , n*
sizeof(
realnum) );
852 memset( &
opac.E2TauAbsOut[lo] , 0 , n*
sizeof(
realnum) );
853 memset( &
opac.TauAbsTotal[lo] , 0 , n*
sizeof(
realnum) );
855 for( i=lo; i <= ihi; i++ )
857 opac.TauTotalGeo[0][i] =
opac.taumin;
859 opac.TauScatGeo[0][i] =
opac.taumin;
861 opac.ExpZone[i] = 1.;
862 opac.E2TauAbsFace[i] = 1.;
863 opac.ExpmTau[i] = 1.;
864 opac.OpacStatic[i] = 1.;
STATIC void fill(double fenlo, double fenhi, double resolv, long int *n0, long int *ipnt, bool lgCount)