/*OpacityCreate compute initial set of opacities for all species */ /****************************************************************************** *NB NB NB NB NB NB NB NB NB NB NB NB NB NB * everything set here must be written to the opacity store files * ****************************************************************************** */ /*lint -e713 loss of precission unsigned int to long */ /*lint -e737 loss of sign in promotion */ #include #include "cddefines.h" #include "taulines.h" #define NCSH2P 10 #include "physconst.h" #include "nhe1lvl.h" #include "opacpoint.h" #include "elmton.h" #include "opac.h" #include "i3727.h" #include "ph1com.h" #include "trace.h" #include "hydrogenic.h" #include "rfield.h" #include "path.h" #include "he.h" #include "nhe1.h" #include "iphmin.h" #include "em2nu.h" #include "nh.h" #include "noptot.h" #include "nxray.h" #include "iphe1l.h" #include "stopcalc.h" #include "phfiton.h" #include "ipoint.h" #include "dustop.h" #include "phfit.h" #include "hypho.h" #include "makeopacity.h" #include "ehe12p.h" #include "hmiopc.h" #include "eva2nu.h" #include "ehe22p.h" #include "hpfit.h" #include "ofit.h" #include "rayleh.h" #include "opacrm.h" #include "opacpl.h" #include "opacity.h" void OpacityCreate(void) { long int i, ipZ, limit, n , nelem; float crs, opart[7]; double dx, eps, thom, thres, x, xnew, z; /* this will be the file name we will try to open for the opacities */ char chFileName[81]; FILE *ioOpac; /* will be pointer to opacity.opc data file */ /* remember whether opacities have ever been evaluated in this coreload */ static int lgOpEval = FALSE; /* fits to cross section for photo dist of H_2^+ */ static float csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f, 6.0f,14.28f,2.7f}; # ifdef DEBUG_FUN fputs( "<+>OpacityCreate()\n", debug_fp ); # endif /* H2+ h2plus h2+ */ /* this is option to stop at certain optical depth */ if( StopCalc.taunu > 0. ) { StopCalc.iptnu = ipoint(StopCalc.taunu); StopCalc.iptnu = MIN2(StopCalc.iptnu,NCELL-1); } else { StopCalc.iptnu = NCELL; StopCalc.tauend = 1e30f; } /* make and print dust opacities * fill in dstab and dstsc, totals, zero if no dust * may be different if different grains are turned on */ dustop(); /* flag lgOpEval says whether opacity stack has been generated * only do this one time per core load */ if( lgOpEval ) { /* this is not the first time code called */ if( trace.lgTrace ) { fprintf( ioQQQ, " OPAC0 called but NOT evaluated.\n" ); } # ifdef DEBUG_FUN fputs( " <->OpacityCreate()\n", debug_fp ); # endif return; } lgOpEval = TRUE; /* should we even try to bring in opacity file? */ if( opac.lgUseFileOpac ) { /* option to compile opacities into file for later use */ /* first try local directory */ strcpy( chFileName , "opacity.opc" ); ioOpac = fopen( chFileName , "rb" ); if( ioOpac == NULL ) { /* did not get it here, lets try the path if it was set */ if( lgDataPathSet ) { /* first stuff the path name itself */ strcpy( chFileName , chDataPath ); strcat( chFileName , "opacity.opc" ); /* now try to open this file */ ioOpac = fopen( chFileName , "rb" ); if( ioOpac == NULL ) { /* still not here, generate new opacities */ opac.lgOpacExist = FALSE; } else { /* found it, will use this file */ opac.lgOpacExist = TRUE; } } else { opac.lgOpacExist = FALSE; } } else { opac.lgOpacExist = TRUE; } /* opac.lgCompileOpac is set TRUE with compile opacities command */ if( !opac.lgCompileOpac && opac.lgOpacExist && (ioOpac!=NULL) ) { /* opacity file exists and have not been asked to compile it */ /* first get the big stack of opacity data */ n = (long)fread( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOpac ); if( ((unsigned)n-sizeof(opac.OpacStack)) != 0 ) { fprintf( ioQQQ, " problem trying to read opacity.opc\n" ); fprintf( ioQQQ, " I expected to read %li words, but fread returned only %li\n", (long)sizeof(opac.OpacStack),n); /* did we hit end of file? */ if( feof(ioOpac ) ) { fprintf( ioQQQ, " end of file hit\n" ); } /* some other error condition? */ if( ferror(ioOpac ) ) { fprintf( ioQQQ, " end of file hit\n" ); } puts( "[Stop in OpacityCreate]" ); exit(1); } /* NG following must exactly mirror contents of common OpacPoint */ /* we will add unwritten bytes to this, then check if still zero at end */ n = 0; n += abs( (long)fread( OpacPoint.ipHOpac, 1,sizeof(OpacPoint.ipHOpac), ioOpac ) -(long)sizeof(OpacPoint.ipHOpac)); n += abs( (long)fread( &OpacPoint.ipRayScat,1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.iopcom, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.ippr, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.ioppr, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.ipBrems, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.iphmra, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.iphmop, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( OpacPoint.ih2pnt, 1,sizeof(OpacPoint.ih2pnt), ioOpac ) - (long)sizeof(OpacPoint.ih2pnt)); n += abs( (long)fread( &OpacPoint.ih2pof, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( OpacPoint.iophe1, 1,sizeof(OpacPoint.iophe1), ioOpac ) - (long)sizeof(OpacPoint.iophe1)); n += abs( (long)fread( &OpacPoint.ioptri, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( OpacPoint.ipElement, 1,sizeof(OpacPoint.ipElement), ioOpac ) - (long)sizeof(OpacPoint.ipElement)); n += abs( (long)fread( OpacPoint.in1, 1,sizeof(OpacPoint.in1), ioOpac ) - (long)sizeof(OpacPoint.in1)); n += abs( (long)fread( OpacPoint.ipo3exc, 1,sizeof(OpacPoint.ipo3exc), ioOpac ) - (long)sizeof(OpacPoint.ipo3exc)); n += abs( (long)fread( OpacPoint.ipo3exc3, 1,sizeof(OpacPoint.ipo3exc3), ioOpac ) - (long)sizeof(OpacPoint.ipo3exc3)); n += abs( (long)fread( OpacPoint.ipo1exc, 1,sizeof(OpacPoint.ipo1exc), ioOpac ) - (long)sizeof(OpacPoint.ipo1exc)); n += abs( (long)fread( &OpacPoint.iopo2d, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.ipmgex, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( &OpacPoint.ipOpMgEx, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); n += abs( (long)fread( OpacPoint.ica2ex, 1,sizeof(OpacPoint.ica2ex), ioOpac ) - (long)sizeof(OpacPoint.ica2ex)); n += abs( (long)fread( &OpacPoint.ica2op, 1,sizeof(long) , ioOpac ) - (long)sizeof(long)); /* following was written as AnuOrg and will be checked against it as a sanity check*/ n += abs( (long)fread( opac.tmn, 1,sizeof(opac.tmn), ioOpac ) - (long)sizeof(opac.tmn)); n += abs( (long)fread( em2nu.Hy2nu, 1,sizeof(em2nu.Hy2nu), ioOpac ) - (long)sizeof(em2nu.Hy2nu)); n += abs( (long)fread( em2nu.he12nu, 1,sizeof(em2nu.he12nu), ioOpac ) - (long)sizeof(em2nu.he12nu)); n += abs( (long)fread( em2nu.he22nu , 1,sizeof(em2nu.he22nu), ioOpac ) - (long)sizeof(em2nu.he22nu)); /* now check whether n is still zero, problem if not */ if( n!= 0 ) { fprintf( ioQQQ, " problem trying to read opacity.opc pointers\n" ); fprintf( ioQQQ, " fread was short by %li\n", n); puts( "[Stop in OpacityCreate]" ); exit(1); } if( opac.tmn[rfield.nupper-1] != rfield.AnuOrg[rfield.nupper-1] ) { fprintf( ioQQQ, " Energy grid in opacity file opacity.opc opacity.pnt do not agree with this version of the code - please recompile.\n" ); fprintf( ioQQQ, " Values are:%12.4e%12.4e\n", rfield.AnuOrg[rfield.nupper-1], opac.tmn[rfield.nupper-1] ); puts( "[Stop in OpacityCreate]" ); exit(1); } if( trace.lgTrace ) { fprintf( ioQQQ, " OPAC0 called, values from file =%s=\n",chFileName ); } /* note that opacities now exist */ opac.lgOpacExist = TRUE; /* zero out opac since this array sometimes addressed before AddOpac called */ for( i=0; i < NCELL; i++ ) { opac.opac[i] = 0.; } # ifdef DEBUG_FUN fputs( " <->OpacityCreate()\n", debug_fp ); # endif return; } else { opac.lgOpacExist = FALSE; } } else { /* this branch if NO FILE OPACITY command entered */ opac.lgOpacExist = FALSE; } if( trace.lgTrace ) { fprintf( ioQQQ, " OPAC0 called, evaluating.\n" ); } /* zero out opac since this array sometimes addressed before AddOpac called */ for( i=0; i <= NCELL-1; i++ ) { opac.opac[i] = 0.; } /* nOpacTot is number of opacity cells OPSV filled so far by OPAC0 */ noptot.nOpacTot = 0; /* hydrogenic photoionization */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { /* this is special pointer to the opacity offset for hydrogenic species only */ OpacPoint.ipHOpac[ipZ][IP1S] = noptot.nOpacTot + 1; /* this is the same number, but in the general purpose pointer array */ OpacPoint.ipElement[2][0][ipZ][ipZ] = OpacPoint.ipHOpac[ipZ][IP1S]; /* the lyman continuum, which extends to the high energy limit of the code */ for( i=nh.ipHn[ipZ][IP1S]-1; i < rfield.nupper; i++ ) { thres = MAX2(rfield.AnuOrg[i]*EVRYD , PH1COM.PH1[0][0][ipZ][0]); phfit(ipZ+1,1,1,thres,&crs,PhFitOn.lgPhFit); opac.OpacStack[i-nh.ipHn[ipZ][IP1S]+OpacPoint.ipHOpac[ipZ][IP1S]] = crs* 1e-18f; } noptot.nOpacTot += rfield.nupper - nh.ipHn[ipZ][IP1S] + 1; /* Balmer, Paschen, etc continua, hbf only goes up to 10, * and NHPLPHOT is defined to 10 in opacpoint.h, * hpfit resolves 2s and 2p */ for( n=IP2S; n < NHPLPHOT; n++ ) { double energy; OpacPoint.ipHOpac[ipZ][n] = noptot.nOpacTot + 1; /* first make sure that first energy point is at least near the limit */ assert( rfield.AnuOrg[nh.ipHn[ipZ][n]-1] > 0.98 * hydro.HLevNIonRyd[ipZ][n] ); for( i=nh.ipHn[ipZ][n]-1; i < nh.ipHn[ipZ][IP1S]; i++ ) { /* for first cell, depending on the curent resolution of the energy mesh, * the center of the first cell can be below the ionization limit of the * level. do not let the energy fall below this limit */ energy = MAX2( hydro.HLevNIonRyd[ipZ][n]*1.001 , rfield.AnuOrg[i] ); opac.OpacStack[i-nh.ipHn[ipZ][n]+OpacPoint.ipHOpac[ipZ][n]] = (float)hpfit(ipZ+1,n,energy*EVRYD); /*(float)hpfit(ipZ+1,n,rfield.AnuOrg[i]*EVRYD);*/ /* make sure it is positive */ assert( opac.OpacStack[i-nh.ipHn[ipZ][n]+OpacPoint.ipHOpac[ipZ][n]] > 0.); } /* this checks that the photo cross sections we just generated agree * with newly computed exact results, at threshold */ # if !defined(NDEBUG) { enum {DEBUG=TRUE}; if( DEBUG && ipZ==0 && n>IP2P ) { float anu[1]={1.} , cs[1] ; double error ; /* photon energy where cross section will be evaluated, * this is in Ryd */ energy = hydro.HLevNIonRyd[ipZ][n]; anu[0] = (float)(energy*1.01); hypho( /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */ (double)(ipZ+1), /* principle quantum number */ n, /* lowest angular momentum */ 0, /* highest angular momentum */ n-1, /* scaled lowest photon energy, * => incorrect?? in units of zed^2/n^2, * at which the cs will be done */ energy, /* number of points to do */ 1, /* array of energies (in units given above) */ anu , /* calculated photoionization cross section in cm^-2 */ cs); error = fabs(cs[0] - opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]-1] )/ ( (cs[0] + opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]-1] )/2.); /*fprintf(ioQQQ,"z=%ld n=%ld error %g old %e new %e\n",ipZ,n, error, opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]-1] ,cs[0] );*/ assert( error < 0.05 ); } } # endif /* possible that center of first cell is below IP of level, * even though high end of level has edge - hpfit returned * zero in this case - update to next higher cell */ if( opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]-1] == 0. ) opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]-1] = opac.OpacStack[OpacPoint.ipHOpac[ipZ][n]]; noptot.nOpacTot += nh.ipHn[ipZ][IP1S] - nh.ipHn[ipZ][n] + 1; if( noptot.nOpacTot > NOPSV - NCELL ) { fprintf( ioQQQ, " hydrogenic opacities need NOPSV larger.\n" ); fprintf( ioQQQ, " ipZ, n, NOPSV=%3ld%3ld%9ld\n", ipZ, n, NOPSV ); puts( "[Stop in OpacityCreate]" ); exit(1); } } } } /* Lyman alpha damping wings - Rayleigh scattering */ OpacPoint.ipRayScat = noptot.nOpacTot + 1; for( i=0; i < nh.ipHn[0][IP1S]; i++ ) { opac.OpacStack[i-1+OpacPoint.ipRayScat] = (float)rayleh(rfield.AnuOrg[i]); } noptot.nOpacTot += nh.ipHn[0][IP1S] - 1 + 1; /* hydrogen two photon emission */ limit = HydroLines[0][IP2P][IP1S].ipCont; xnew = eva2nu(rfield.AnuOrg[0]); em2nu.Hy2nu[0] = (float)(xnew*rfield.widflx[0]); for( i=1; i < limit; i++ ) { xnew = eva2nu(rfield.AnuOrg[i]); em2nu.Hy2nu[i] = (float)(xnew*rfield.widflx[i]); } /* HeI two photon emission */ for( i=0; i < iphe1lCom.iphe1l[0][1]; i++ ) { em2nu.he12nu[i] = (float)(ehe12p(rfield.AnuOrg[i])*rfield.widflx[i]); } /* HeII two photon emission */ limit = HydroLines[1][IP2P][IP1S].ipCont; for( i=0; i < limit; i++ ) { em2nu.he22nu[i] = (float)(ehe22p(rfield.AnuOrg[i])*rfield.widflx[i]); } /* assume Thomson scattering up to ipCKshell, 20.6 Ryd=0.3 keV */ thom = 6.65e-25; OpacPoint.iopcom = noptot.nOpacTot + 1; for( i=0; i < nxray.ipCKshell; i++ ) { opac.OpacStack[i-1+OpacPoint.iopcom] = (float)thom; } /* Klein-Nishina from eqn 7.5, Rybicki and Lightman */ for( i=nxray.ipCKshell; i < rfield.nupper; i++ ) { dx = rfield.AnuOrg[i]/3.7573e4; opac.OpacStack[i-1+OpacPoint.iopcom] = (float)(thom*3.e0/4.e0*((1.e0 + dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+ 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0* dx)/POW3(1.e0 + 2.e0*dx))); } noptot.nOpacTot += rfield.nupper - 1 + 1; /* pair production */ OpacPoint.ioppr = noptot.nOpacTot + 1; for( i=OpacPoint.ippr-1; i < rfield.nupper; i++ ) { /* pair production heating rate for unscreened H + He * fit to figure 41 of Experimental Nuclear Physics, * Vol1, E.Segre, ed */ x = rfield.AnuOrg[i]/7.512e4*2.; opac.OpacStack[i-OpacPoint.ippr+OpacPoint.ioppr] = (float)(5.793e-28* POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. + x*(0.130471301 + x*0.000524906)))); } noptot.nOpacTot += rfield.nupper - OpacPoint.ippr + 1; /* brems (free-free) opacity */ OpacPoint.ipBrems = noptot.nOpacTot + 1; for( i=0; i < rfield.nupper; i++ ) { /* missing factor of 1E-20 to avoid underflow * free free opacity needs g(ff)*(1-exp(hn/kT))/SQRT(T)*1E-20 */ opac.OpacStack[i-1+OpacPoint.ipBrems] = (float)(1.03680e-18/POW3(rfield.AnuOrg[i])); } noptot.nOpacTot += rfield.nupper - 1 + 1; OpacPoint.iphmra = noptot.nOpacTot + 1; for( i=0; i < rfield.nupper; i++ ) { /* following is ratio of h minus to neut h opacity */ opac.OpacStack[i-1+OpacPoint.iphmra] = (float)(0.1175*rfield.anusqr[i]); } noptot.nOpacTot += rfield.nupper - 1 + 1; OpacPoint.iphmop = noptot.nOpacTot + 1; for( i=iphminCom.iphmin-1; i < nhe1Com.nhe1[0]; i++ ) { /* H- hminus H minus bound-free opacity */ opac.OpacStack[i-iphminCom.iphmin+OpacPoint.iphmop] = (float)(hmiopc(rfield.AnuOrg[i])); } noptot.nOpacTot += nhe1Com.nhe1[0] - iphminCom.iphmin + 1; /* H2+ H2P h2plus photoabsorption * cs from Buckingham Reid and Spence, MN 112m 382, 0 K temp */ opacrm(OpacPoint.ih2pnt[0],OpacPoint.ih2pnt[1],csh2p,NCSH2P,&OpacPoint.ih2pof, "H2+ "); /* Helium, triplet 23S from Steward 1979, quoted in Mathisen compendium */ opacpl(he.nhei3,nhe1Com.nhe1[0],5.5e-18,1.35,&OpacPoint.ioptri); /* HeI singlets neutral helium ground */ OpacPoint.iophe1[0] = noptot.nOpacTot + 1; OpacPoint.ipElement[2][0][0][1] = OpacPoint.iophe1[0]; for( i=nhe1Com.nhe1[0]-1; i < rfield.nupper; i++ ) { phfit(2,2,1,rfield.AnuOrg[i]*EVRYD,&crs,PhFitOn.lgPhFit); opac.OpacStack[i-nhe1Com.nhe1[0]+OpacPoint.iophe1[0]] = (float)(crs*1e-18); } noptot.nOpacTot += rfield.nupper - nhe1Com.nhe1[0] + 1; /* HeI neutral Helium 21s from Stewart JPhysB 11, L431 */ opacpl(nhe1Com.nhe1[1],nhe1Com.nhe1[0],0.4*8.7e-18,1.5,&OpacPoint.iophe1[1]); /* HeI helium singlet bound free excited states */ z = 1.; for( n=3; n <= NHE1LVL; n++ ) { crs = (float)(7.906e-18*(float)(n)/(z*z)); opacpl(nhe1Com.nhe1[n-1],nhe1Com.nhe1[0],crs,3.0,&OpacPoint.iophe1[n-1]); } /* these are opacity offset points that would be defined in MakeOpacity, * but this routine will not be called for H and He * generate all heavy element opacities, everything heavier than He, * nelem is on the C scale, so Li is 2 */ /*>>chng 99 jan 27, do not reevaluate hydrogenic opacity below */ for( nelem=2; nelem < LIMELM; nelem++ ) { MakeOpacity(nelem); } /* now add on some special cases, where exicted states, etc, come in */ /* Nitrogen * >>refer Henry, R., ApJ 161, 1153. * photoionization of excited level of N+ */ opacpl(OpacPoint.in1[0],OpacPoint.in1[1],9e-18,1.75,&OpacPoint.in1[2]); /* atomic Oxygen * only do this if 1996 Verner results are used */ if( !PhFitOn.lgPhFit ) { /* integrate over energy range of the valence shell of atomic oxygen*/ for( i=OpacPoint.ipElement[0][2][0][7]-1; i < OpacPoint.ipElement[1][2][0][7]; i++ ) { /* call special routine to evaluate partial cross section for OI shells */ eps = rfield.AnuOrg[i]*EVRYD; ofit(eps,&crs,opart); /* this will be total cs of all processes leaving shell 3 */ crs = opart[0]; for( n=1; n < 6; n++ ) { /* add up table of cross sections */ crs += opart[n]; } /* convert to cgs and overwrite cross sections set by MakeOpacity */ crs *= 1e-18f; opac.OpacStack[i-OpacPoint.ipElement[0][2][0][7]+OpacPoint.ipElement[2][2][0][7]] = crs; } } /* Henry nubmers for 1S excit state of OI, OP data very sparse */ opacpl(OpacPoint.ipo1exc[0],OpacPoint.ipo1exc[1],4.64e-18,0.,&OpacPoint.ipo1exc[2]); /* photoionization of excited level of O2+ 1D making 5007 * fit to TopBase Opacity Project cs */ opacpl(OpacPoint.ipo3exc[0],OpacPoint.ipo3exc[1],3.8e-18,0.,&OpacPoint.ipo3exc[2]); /* photoionization of excited level of O2+ 1S making 4363 */ opacpl(OpacPoint.ipo3exc3[0],OpacPoint.ipo3exc3[1],5.5e-18,0.01, &OpacPoint.ipo3exc3[2]); /* photoionization to excited states of O+ */ OpacPoint.iopo2d = noptot.nOpacTot + 1; thres = rfield.AnuOrg[i3727.i2d-1]; for( i=i3727.i2d-1; i < nh.ipHn[1][IP1S]; i++ ) { crs = (float)(3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38* pow(rfield.AnuOrg[i]/thres,-2.5))); opac.OpacStack[i-i3727.i2d+OpacPoint.iopo2d] = crs; } noptot.nOpacTot += nh.ipHn[1][IP1S] - i3727.i2d + 1; /* magnesium * photoionization of excited level of Mg+ * fit to opacity project data Dima got */ OpacPoint.ipOpMgEx = noptot.nOpacTot + 1; for( i=OpacPoint.ipmgex-1; i < nh.ipHn[0][IP1S]; i++ ) { opac.OpacStack[i-OpacPoint.ipmgex+OpacPoint.ipOpMgEx] = (float)((0.2602325880970085 + 445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))* 1e-18); } noptot.nOpacTot += nh.ipHn[0][IP1S] - OpacPoint.ipmgex + 1; /* Calcium * excited states of Ca+ */ opacpl(OpacPoint.ica2ex[0],OpacPoint.ica2ex[1],4e-18,1.,&OpacPoint.ica2op); if( trace.lgTrace ) { fprintf( ioQQQ, " OPAC0 return OK, number of opacity cells used in OPSC=%7ld and OPSV is dimensioned%7ld\n", noptot.nOpacTot, NOPSV ); } if( NOPSV - noptot.nOpacTot < NCELL ) { fprintf( ioQQQ, " Warning, only%4ld cells left in OPSV. Increase NOPSV.\n", NOPSV - noptot.nOpacTot ); } /* option to compile opacities into file for later use * this is executed if the 'compile opacities' command is entered */ if( opac.lgCompileOpac ) { ioOpac = fopen("opacity.opc","wb"); if( ioOpac == NULL ) { fprintf( ioQQQ, " problem trying to open opacity.opc\n" ); puts( "[Stop in OpacityCreate]" ); exit(1); } /* write the data, remembering how many bytes were dumped */ n = fwrite( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOpac ); if( (n-sizeof(opac.OpacStack)) != 0 ) { fprintf( ioQQQ, " problem trying to write opacity.opc\n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite returned only %li\n", (long)sizeof(opac.OpacStack),n); puts( "[Stop in OpacityCreate]" ); exit(1); } /* NG following must exactly mirror contents of common OpacPoint */ /* we will add unwritte bytes to this, then check if still zero at end */ n = 0; n += abs( (int)fwrite( OpacPoint.ipHOpac, 1,sizeof(OpacPoint.ipHOpac), ioOpac ) -sizeof(OpacPoint.ipHOpac)); n += abs( fwrite( &OpacPoint.ipRayScat,1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.iopcom, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.ippr, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.ioppr, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.ipBrems, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.iphmra, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.iphmop, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( OpacPoint.ih2pnt, 1,sizeof(OpacPoint.ih2pnt), ioOpac ) - sizeof(OpacPoint.ih2pnt)); n += abs( fwrite( &OpacPoint.ih2pof, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( OpacPoint.iophe1, 1,sizeof(OpacPoint.iophe1), ioOpac ) - sizeof(OpacPoint.iophe1)); n += abs( fwrite( &OpacPoint.ioptri, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( OpacPoint.ipElement, 1,sizeof(OpacPoint.ipElement), ioOpac ) - sizeof(OpacPoint.ipElement)); n += abs( fwrite( OpacPoint.in1, 1,sizeof(OpacPoint.in1), ioOpac ) - sizeof(OpacPoint.in1)); n += abs( fwrite( OpacPoint.ipo3exc, 1,sizeof(OpacPoint.ipo3exc), ioOpac ) - sizeof(OpacPoint.ipo3exc)); n += abs( fwrite( OpacPoint.ipo3exc3, 1,sizeof(OpacPoint.ipo3exc3), ioOpac ) - sizeof(OpacPoint.ipo3exc3)); n += abs( fwrite( OpacPoint.ipo1exc, 1,sizeof(OpacPoint.ipo1exc), ioOpac ) - sizeof(OpacPoint.ipo1exc)); n += abs( fwrite( &OpacPoint.iopo2d, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.ipmgex, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( &OpacPoint.ipOpMgEx, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( OpacPoint.ica2ex, 1,sizeof(OpacPoint.ica2ex), ioOpac ) - sizeof(OpacPoint.ica2ex)); n += abs( fwrite( &OpacPoint.ica2op, 1,sizeof(long) , ioOpac ) - sizeof(long)); n += abs( fwrite( rfield.AnuOrg, 1,sizeof(rfield.AnuOrg), ioOpac ) - sizeof(rfield.AnuOrg)); n += abs( fwrite( em2nu.Hy2nu, 1,sizeof(em2nu.Hy2nu), ioOpac ) - sizeof(em2nu.Hy2nu)); n += abs( fwrite( em2nu.he12nu, 1,sizeof(em2nu.he12nu), ioOpac ) - sizeof(em2nu.he12nu)); n += abs( fwrite( em2nu.he22nu , 1,sizeof(em2nu.he22nu), ioOpac ) - sizeof(em2nu.he22nu)); /* now check whether n is still zero, problem if not */ if( n!= 0 ) { fprintf( ioQQQ, " problem trying to write opacity.opc pointers\n" ); fprintf( ioQQQ, " fwrite was short by %li\n", n); puts( "[Stop in OpacityCreate]" ); exit(1); } fclose(ioOpac); fprintf( ioQQQ, " Compile opacity completed ok - stopping.\n" ); puts( "[Stop in OpacityCreate]" ); exit(1); } # ifdef DEBUG_FUN fputs( " <->OpacityCreate()\n", debug_fp ); # endif return; } /*lint +e713 loss of precission unsigned int to long */ /*lint +e737 loss of sign in promotion */