/* smpopac2.c, Apr 2006:  This is like simpopac.c except for two changes:
  (1) The set of output wvl or photon-energy intervals has been changed. 
    Now the output list extends beyond the Lyman limit at _nu_ = 10.97. 
  (2) The input-data interval in _u_ has been reduced by a factor of
    two (i.e., finer spacing) with pseudo-cubic interpolation. 
    This requires two large storage arrays nu[] and sss[], allocated
    by malloc(). 	
   --------------------------------------------------------------------- */

/* simpopac.c,  April 2006:  produces ascii list of simplified opacities 
   kappa(nu), using input files produced by the Opacity Project "monop.out" 
   program.  (This program won't work if some details of the input-file 
   format are altered.)  We assume that opacities beyond the Lyman limit 
   are very high (not calculated).   
   Output quantity _nu_ is (1 micron)/wvl.
   ----------------------------------------------------------------- */

#include <stdlib.h>  
#include <stdio.h>  
#include <math.h> 
#include <fcntl.h>

main()
 {
  char str[400],fnameout[31],fnamein[11][31],datastr[4200];
  int i,j,k,l,m,n,o,p,place,n_ionstates,offset;  
  int fdout,fdin,ite[11],jne,izz[11],datastrsize;
  int nelems,elemn,ninpts[11];   
  int nintervals,intrvln,m1[12],dm12[12];  
  double a,b,c,d,e,f,g,h,nux;
  double atom_mass[11],rosseland[11],fion[11][10],abund[11]; 
  double e_per_n[11];  
  double kappa_sc,kappa_rm,sum_rm,kappa_nu;    
  double log10temp,temp,log10ne,dens_atom[12],dens_e,rho;
  double *nu,*sss,sss_min;    
  double nulim[135],nu_c[134],dnu[134],u,du,ufac;   


  
    /* Prepare the big storage arrays so that values other than 
       those to be read in will be clearly out of bounds: */
  nu  = malloc(2112000);  /* (11 elements)*(24000 samples)*(8 bytes) */ 
  sss = malloc(2112000);   
  for (j=0;j<264000;j++) { nu[j] = 100.;  sss[j] = 1.e20; }

  printf("\n\n\n-------------------------------------------------------\n");
  fdout = -1;
  while (fdout < 0)
   {
    printf("\nOutput file name (new ascii file) :  "); 
    gets(str);  if (str[0] == 'q' && str[1] == '\0') exit(0);
    fnameout[30] = '\0';  
    for (i=0;i<30;i++)
      { fnameout[i] = str[i];  if (fnameout[i] == '\0') break; } 
    j = open(fnameout,0);
    if (j >= 0)
      { printf("\n***  \"%s\" already exists  ***\n",fnameout);  close(j); } 
    else
     { 
      fdout= open(fnameout,O_CREAT|O_EXCL|O_WRONLY,0666);
      if (fdout < 0) 
         printf("\n***  Can\'t create \"%s\"  ***\n",fnameout); 
     }
   }
  printf("Output file is \"%s\".\n",fnameout);  
  sprintf(str,"---\nFile produced by simpopac.c:  "); 
  for (i=0;i<100;i++) if (str[i] == '\0') break;
  write(fdout,str,i);  

  k = 0;  
  while (k != 1)
   {
    printf("\nInteger jne = 4*log10(ne) :  ");
    gets(str); 
    if (str[0] == 'q') 
     { 
      printf("\n---  Remember to delete \"%s\"  ***\n",fnameout);
      close(fdout);  exit(0); 
     }
    k = sscanf(str,"%d",&jne);   
    if (k == 1)
     {
      log10ne = 0.25*((double) jne);   
      while (str[0] != 'y' && str[0] != 'q' && str[0] != 'n')  
       {
        printf("jne = %d,  log10(ne) = %.4lf;  OK?  (y/n) :  ",jne,log10ne);  
        gets(str);
       }
      if (str[0] == 'q') 
       { 
        printf("\n---  Remember to delete \"%s\"  ***\n",fnameout);
        close(fdout);  exit(0); 
       }
      if (str[0] != 'y') k = 0; 
     }
   }

  dens_e = exp(2.302585093*log10ne);  /* electron density */    
  for (elemn=0;elemn<11;elemn++) e_per_n[elemn] = 0.;  


  /*  Define 134 intervals in _nu_:  */
  nulim[0] = 0.17;    
  for (j=1;j<=18;j++)   nulim[j] = 1.1*nulim[j-1];     /* 0.187--0.945 */ 
  for (j=19;j<=44;j++)  nulim[j] = 1.05*nulim[j-1];    /* 0.992--3.361 */
  for (j=45;j<=97;j++)  nulim[j] = nulim[j-1] + 0.17;  /* 3.531--12.37 */     
  for (j=98;j<=134;j++) nulim[j] = nulim[j-1] + 0.25;  /* 12.62--21.87 */
  for (j=0;j<=133;j++)   
   { 
    nu_c[j] = 0.5*(nulim[j] + nulim[j+1]);  
    dnu[j]  = nulim[j+1] - nulim[j];  
   }

  printf("\nNote:  When this program queries for the relative abundance of ");
  printf("each\nelement, it means by number (not by mass) but the ");
  printf("normalization can\nbe arbitrary.  (The program will take care ");
  printf("of normalization.)\nMax Nelements = 11.\n"); 

  nelems = 0;  
  for (elemn=0;elemn<11;elemn++)   /* This is a long loop */  
   {

    datastrsize = 0;   
    fdin = -1; 
       /* Preliminary information in the input file for _elemn_ 
          is entirely read in the following "while fdin < 0" 
          procedure: */
    while (fdin < 0)   
     {

      printf("\nInput file %d (or \'x\' if finished) :  ",elemn+1);
      gets(str);
      if (str[0] == 'q' && str[1] == '\0')
       { 
        printf("\n--- Remember to delete new file \"%s\" ---\n",fnameout);
        close(fdout);  exit(0);
       }
      if (str[0] == 'x' && str[1] == '\0' && elemn) break; 
           /* (above) break from "while fdin < 0" */  

      fnamein[elemn][30] = '\0';  
      for (p=0;p<30;p++) 
        { fnamein[elemn][p] = str[p];  if (str[p] == '\0') break; }

      fdin = open(&fnamein[elemn][0],0);    

      if (fdin < 0) 
        printf("\n***  Can\'t open \"%s\"  ***\n",&fnamein[elemn][0]);  
      else   
       { 
        datastrsize = read(fdin,datastr,4000);  
            if (datastrsize < 0) datastrsize = 0;
        if (datastrsize < 100)
         { 
          printf("\n***  %d bytes, too small  ***\n",datastrsize);  
          datastrsize = 0;  
         }
       }

           /* Required structure of input file for each chem element:
                  (file produced by Opacity Project program 
                   "monop.out".  Each text line is terminated by 
                   a LF character, byte value = decimal 10.) 
              1st text line:  Z, ite, atommass, aseveral other quantities;
              2nd line:  several parameters, not used here;
              3rd line:  two integers i j identifying significant
                ionization states, see below.
              Lines 4 ... 4+j-i  (i.e., j-i+1 lines):  On each line,
                at integer k and an ionization fraction.  The integer
                identifies the ionization stage in a peculiar way,
                i.e.  net charge z = Z - k - 1.   (z = -1, the negative
                ion, can be important for hydrogen.)      
              So far there are j-i+4 lines of preliminary data, where
                i and j are listed on line 3.  
              Every line after this lists two f.p. numbers, _nu_ and
                _s_.    The _nu_ parameter is h*nu/k*T.  Quantity _s_
                is  (cross-section)/{1 - exp(-u)} where "cross-section"
                means the cross-section per atom of this element 
                (including all ionization states).  Contributions of
                this element to the Thomson-scattering and free-free
                opacties are included.  Note that the correction for
                stimulated emission, a factor of {1 - exp(-u)}, has
                NOT been made in the listed values of _s_;  in other
                words, to get effective cross-sections we must reduce
                each value by that factor.  Moreover, the listed
                cross-sections are in units of a0^2 (Bohr radius
                squared), not cm^2.   The Thomson scattering cross-
                section for an electron is about 2.376e-8*a0^2.
                Unfortunately the cross-section data are listed to
                a precision of only three decimal places, we must
                be careful about rounding errors in cases where
                Thomson scattering accounts for almost all of the
                opacity (otherwise we could get illusory negative
                values of kappa_abs).

              Throughout the following procedure, _datastrsize_ is the
              number of useful ascii bytes currently in str[].  Quantity
              _place_ is a current location in str[] immediately following
              an end-of-line character (LF).  Unfortunately the procedure
              is cumbersome because we must read through 6 or more text
              lines of prelminary information before we get to the main
              data list (see format description above).   

              If any obvious errors or inconsistencies are found,
              datastrsize will be reset to zero, which will force
              the "while fdin < 0" procedure to be repeated. 
              ----------------------------------------------------  */ 


      if (datastrsize)  /* Read the first text line */        
       {
        k = sscanf(datastr,"%d%d%le",&izz[elemn],&ite[elemn],&atom_mass[elemn]);  
        if (k == 3)
         {
          for (j=0;j<elemn;j++)   /* compare Z with previous entries */ 
           {
            if (izz[j] == izz[elemn]) 
             {
              printf("\n***  izz = %d duplicates a preceding input file  ***\n",izz[elemn]);  
              datastrsize = 0;
             }
           } 
         }
        else  /* i.e., if k != 3 */
         { 
          printf("\n***  Can\'t read first line  ***\n"); 
          datastrsize = 0; 
         }  
       }  

        /* Temperature parameter _ite_ is specified on 1st line: */
      if (datastrsize)   
       {
        if (!elemn)   /* first specified element, get temperature */
         { 
          log10temp = 0.025*((double) ite[0]); 
          temp = exp(2.302585093*log10temp);  
          ufac = 14388.35/temp;       
            /* (above)  u = h*nu/k*T = nu/ufac;  
               _nu_ is (1 micron)/wvl  in this program.  */ 
          sprintf(str,"\nlog10(T) = %.3lf, log10(ne) = %.1lf  (ite=%d, jne=%d);\n",log10temp,log10ne,ite[0],jne);  
          for (i=0;i<400;i++) if (str[i] == '\0') break;
          write(fdout,str,i);
         } 
        else     /* if elemn > 0, verify that temperature is consistent */  
         {
          if (ite[elemn] != ite[0])   
           { 
            printf("\n***  Wrong ite, %d, should be %d  ***\n",ite[elemn],ite[0]);  
            datastrsize = 0; 
           } 
         }
       }

      if (datastrsize)   /* Read the second line of the input file */ 
       { 
        for (place=0;place<400;place++) 
            if (datastr[place] == 10) break;  /* first LF */ 
        ++place; 
        k = sscanf(&datastr[place],"%le%le%le",&a,&b,&rosseland[elemn]);  
        if (k != 3) 
         { 
          printf("\n***  2nd line doesn't have 3 quantities  ***\n");
          datastrsize = 0; 
         }  
       } 

      if (datastrsize)   /* Read third input line */  
       {
        for (m=place;m<datastrsize;m++) 
           { if (datastr[m] == 10) break; } /* 2nd LF */ 
        place = m + 1;   
        k = sscanf(&datastr[place],"%d%d",&i,&j);   
        n_ionstates = 0;
        if (k == 2)  n_ionstates = j - i + 1;  
        if (k != 2 || n_ionstates < 1 || n_ionstates > 12)  
         {
          printf("\n***  Trouble on 3rd line  ***\n");
          datastrsize = 0; 
         } 
       }

      if (datastrsize)
       {
          /* fion[elemn][i] will be fraction of that atom with
             net charge z = i - 1.  The special case i = 0
             allows for negative hydrogen ions, z = -1. */  

        for (i=0;i<10;i++) fion[elemn][i] = 0.; 

        for (o=0;o<n_ionstates;o++)
         {
          for (m=place;m<datastrsize;m++) { if (datastr[m] == 10) break; }  
          place = m + 1;   /* LF number o+3 */  
          k = sscanf(&datastr[place],"%d%le",&j,&a);   
          if (k == 2)
           { 
            i = izz[elemn] - j;   
            if (i >= 0 && i < 10)  fion[elemn][i] = a; 
            if (a < 0. || a > (1. + 1.e-13)) k = 0; 
           } 
          if (k != 2)
           { 
            printf("\n***  Trouble on line %d  ***\n",o+4);   
            datastrsize = 0;  break;  
           }
         }
       }

      if (datastrsize)  /* renormalize fion[][] */ 
       { 
        a = fion[elemn][0];  for (i=1;i<10;i++) a = a + fion[elemn][i];  
        if (a > 0.995 && a < 1.005)
         { for (i=0;i<10;i++) fion[elemn][i] = fion[elemn][i]/a;  }
        else 
         { 
          printf("\n***  Warning:  Sum of ionization fractions = %.4le,",a);  
          printf("\n     they haven't been renormalized so sum = 1.");
          printf("\n     Hit <Enter> to continue :  ");
          gets(str);   
          if (str[0] == 'q')
           {
            printf("\n---  Remember to delete \"%s\"  ---\n",fnameout);
            close(fdout);  exit(0);
           }
         }
        e_per_n[elemn] = -fion[elemn][0]; 
        for (i=2;i<10;i++)
          e_per_n[elemn] = e_per_n[elemn] + fion[elemn][i]*((double) (i - 1));  
        printf("%.6lf electrons per nucleus of this element.",e_per_n[elemn]);         
       } 
    
        /* At this point, the preliminary information in the input
           file should have been read.  If anything is obviously
           wrong, datastrsize has been reset to zero.  In that 
           case, the following line will force the "while fdin < 0"
           to be repeated, i.e., the input file will be requested
           again. */  
      if (!datastrsize) { close(fdin); fdin = -1; } 

     }   /* --- end of "while fdin < 0" --- */

    if (fdin < 0) break;    
       /* (above) Break from elemn-loop.  This occurs if 'x' was 
          entered instead of an input file name, indicating that 
          the list of elements is complete. */   

    /* Specify the abundance of this element:  (arbitrary relative
       units, abundance by number rather than by mass) */
    k = 0;
    while (k != 1)
     {
      printf("\nRelative abundance of Z = %d :  ",izz[elemn]);
      gets(str);
      if (str[0] == 'q')
       {
        printf("\n---  Remember to delete \"%s\"  ---\n",fnameout);
        close(fdout);    close(fdin);  exit(0);  
       } 
      k = sscanf(str,"%le",&abund[elemn]);
      if (abund[elemn] < 0.) k = 0; 
     }

    sprintf(str,"---\nInput file \"%s\":\niz = %d, m = %.4lf,  rel abund = %.5le.\n",&fnamein[elemn][0],izz[elemn],atom_mass[elemn],abund[elemn]); 
    for (i=0;i<200;i++) if (str[i] == '\0') break;
    write(fdout,str,i);

    sprintf(str,"Ionization (0 is the neutral atom):\n"); 
    for (i=0;i<200;i++) if (str[i] == '\0') break;
    write(fdout,str,i);
    for (o=0;o<10;o++) if (fion[elemn][o]) break; 
    for (p=9;p>=0;p--) if (fion[elemn][p]) break;
    for (j=o;j<=p;j++) 
     { 
      sprintf(str,"%5d  %.7le\n",j-1,fion[elemn][j]);
      for (i=0;i<200;i++) if (str[i] == '\0') break;
      write(fdout,str,i);
     }
    sprintf(str,"%.6lf electrons per nucleus of this element.  \n",e_per_n[elemn]);
    for (i=0;i<200;i++) if (str[i] == '\0') break;
    write(fdout,str,i);

	sss_min = 6.653e-25*e_per_n[elemn];   
    offset = 24000*elemn;    
	ninpts[elemn] = 0;  

         /* copy data into arrays nu[] and sss[]: */

    for (ninpts[elemn]=0;ninpts[elemn]<11994;++ninpts[elemn])   
     {
         /* If necessary, read new data from the current input file: */
        p = datastrsize - place;
        if (p < 60)   
         {
          for (i=0;i<p;i++)  datastr[i] = datastr[place+i]; 
          place = 0;  
          o = read(fdin,&datastr[p],4000);  if (o < 0) o = 0; 
          datastrsize = p + o;  
            /* (Next) Trick to force an end to the data reading when 
               the last input line has been reached:  We do this by 
               making some  text lines that have no numbers.  */ 
          if (o != 4000)  /* end of this input file */
           {
            close(fdin);  fdin = -1;   
            for (i=0;i<100;i++) datastr[datastrsize+i] = ' ';  
            for (i=4;i<100;i=i+5) datastr[datastrsize+i] = 10; 
            datastrsize = datastrsize + 100;
           }
         }

      for (m=place;m<datastrsize;m++) 
      if (datastr[m] == 10) break;  /* find next 10 = LF */
      place = m + 1; 

      k = sscanf(&datastr[place],"%le%le",&u,&b);
      if (k == 2) { if(u <= 0. || b <= 0.) k = 0; }

      if (k != 2)  break;  
	     /* (Above) Have reached end of available data for 
		     this element, break from the ninpts-loop. */ 

      n = offset + ninpts[elemn];  /* address in big array */  
      nu[n] = u/ufac;     /* _u_ is h*nu/k*T  */
      sss[n] = 2.8052e-17*b*(1. - exp(-u));  
          /* (above) The numerical value is a0^2, the square of
             the Bohr radius;  the other factor is the correction    
             for stimulated emission.  Quantity sss[][] is the
             effective cross section per atom of type _elemn, 
             measured in cm^2.   
             There is a small numerical discrepancy.  Strictly
             speaking the numerical constant should be 
             2.8003e-17 cm^2 instead of 2.8052e-17;  but
             that led to some opacities smaller than the
             Thomson scattering value.  The discrepancy is
             less than 0.2 percent, and may result from
             rounding in the OPAC tables.  */  

      if (ninpts[elemn] > 0)
       {
        if (nu[n] <= nu[n-1])   /* must be definitely monotonic */
		  { k = 0;   nu[n] = 100.;  sss[n] = 1.e20;  break; }  
       }

     }  /* ninpts-loop, which copies data to arrays nu[][] and sss[][] */

    if (ninpts[elemn] < 3) 
     {
      printf("\n***  Fewer than 3 input points  ***\n");
      sprintf(str," ^^^ ERR, ignore this element; \n---\n");  
      for (i=0;i<200;i++) if (str[i] == '\0') break;  
      write(fdout,str,i);
      --elemn;
     }
    else  /* Normal case, at least 3 input points */  
     {
      nelems = elemn + 1;  

      if (ninpts[elemn] > 11990) 
       {
        sprintf(str," ^^^ WARNING:  This data list may be incomplete.\n---\n");
        for (i=0;i<200;i++) if (str[i] == '\0') break;
        write(fdout,str,i);  
        printf("\n***  WARNING:  Input list is too long, "); 
        printf("\n     data list may be truncated. ");
        printf("\n     <Continue> :  ");
        gets(str);  
        if (str[0] == 'q') { close(fdout); close(fdin); exit(0); } 
        break;
       } 

        /* Expand input list by interpolation.  First move
		   points {0,1,...ninpts-1} to {1,3,...2*ninpts-1}:  */   
      for (n=ninpts[elemn]-1;n>=0;n--)     
       {
        o = offset + n;  p = offset + 2*n + 1;       
        nu[p] = nu[o];  sss[p] = sss[o];  
	   }   
        /* Redefine total number of stored points that we'll
		   have after extension and interpolation: */ 
      ninpts[elemn] = 2*ninpts[elemn] + 1;
        /* Relative to the offset, the data points copied 
		   from the input file are now in locations 
		   1,3,...ninpts-2.   Create new endpoints at 0
	       and ninpts-1, which definitely enclose the 	   
           entire range spanned by output list nulim[]: */
      a = 0.95*nu[offset+1];   
      if (a >= nulim[0]) a = 0.95*nulim[0];  
      nu[offset] = a;  sss[offset] = 1.05*sss[offset+1];  
      n = offset + ninpts[elemn] - 1;  
      b = nu[n-1] + 0.01; 
      if (b <= nulim[135]) b = nulim[135] + 0.01;  
      nu[n] = b;  sss[n] = sss[n-1];    
        /* Simple interpolation for pts 2 and ninpts-3: */
      n = offset + 2;
      nu[n] = sqrt(nu[n-1]*nu[n+1]);
      sss[n] = sqrt(sss[n-1]*sss[n+1]);   
      n = offset + ninpts[elemn] - 3;  
      nu[n] = 0.5*(nu[n-1] + nu[n+1]);
      sss[n] = sqrt(sss[n-1]*sss[n+1]);
        /* Fancier interpolation for points
		   4, 6, ... ninputs-5: */   
      for (m=4;m<=ninpts[elemn]-5;m=m+2)  
       {
        n = offset + m;  
        nu[n] = 0.5*(nu[n-1] + nu[n+1]);  
        h = (nu[n-3] - nu[n])/(nu[n+1] - nu[n]);   
        g = 0.5*h/(h - 1.);  f = 0.5*h/(h + 1.);   
        e = 1. - f - g;  
        a = e*log(sss[n-3]) + f*log(sss[n-1]) + g*log(sss[n+1]);   
        h = (nu[n+3] - nu[n])/(nu[n+1] - nu[n]); 
        e = 0.5*h/(h + 1.);  f = 0.5*h/(h - 1.);
        g = 1. - e - f; 
        b = e*log(sss[n-1]) + f*log(sss[n+1]) + g*log(sss[n+3]);  
        sss[n] = exp(0.5*(a + b));     
       }

      /* At any wavelength, the opacity must not be less than
	     the amount due to Thomson scattering by free electrons 
		 provided by this element: */   
      for (m=0;m<ninpts[elemn];m++)  
       {
        n = offset + m;  
        if (sss[n] < sss_min) sss[n] = sss_min;  
       }

     }    /* --- end of "else" (ninpts >= 3) --- */
      
   }     /* ---  elemn-loop  --- */


                 /* temporary */
                 printf("\nNelems = %d;   ",nelems);  gets(str); 

      /* At this point, the relevant input data have been stored 
         in arrays nu[] and sss[].   (The nu-values have
         been calculated from the input values of U = h*nu/k*T.) */

  /* ===-=== Next, estimate the atom (nucleus) densities, then save 
     some information in the output file.  Begin by figuring out
     the normalization factor for the specified abundance values,
	 which implies the mass density:  */   
  for (elemn=0;elemn<11;elemn++) dens_atom[elemn] = 0.;  
  a = 0.;  /* Count free electrons per atom abundance */
  for (elemn=0;elemn<nelems;elemn++) a = a + abund[elemn]*e_per_n[elemn];  
  b = dens_e/a;   rho = 0.; 
  for (elemn=0;elemn<nelems;elemn++) 
   {
    dens_atom[elemn] = b*abund[elemn];  
    rho = rho + atom_mass[elemn]*dens_atom[elemn];   
   }
  rho = 1.66053e-24*rho;              /* amu */  
  kappa_sc = 6.6525e-25*dens_e/rho;   /* Thomson scattering */  

  c = log10(rho);
  printf("\nlog10(rho) = %.5lf;\n",c); 

  sprintf(str,"---\nrho = %.6le,     log10(rho) = %.5lf.\n",rho,c);
  for (i=0;i<200;i++) if (str[i] == '\0') break; 
  write(fdout,str,i); 

  sprintf(str,"  n_e    = %.5le,   kappa_sc = %.6lf;\n",dens_e,kappa_sc);
  for (i=0;i<200;i++) if (str[i] == '\0') break; 
  write(fdout,str,i); 
  for (elemn=0;elemn<nelems;elemn++)
   { 
    sprintf(str,"n(iz=%2d) = %.5le  \n",izz[elemn],dens_atom[elemn]);
    for (i=0;i<200;i++) if (str[i] == '\0') break; 
    write(fdout,str,i); 
   }
  sprintf(str,"---\n");
  for (i=0;i<200;i++) if (str[i] == '\0') break; 
  write(fdout,str,i); 

  sprintf(str,"    nu        kappa     k_abs/k_sc  \n");
  for (i=0;i<200;i++) if (str[i] == '\0') break; 
  write(fdout,str,i); 



  /* ===-===  Do the main calculations: */ 
  for (intrvln=0;intrvln<nintervals;intrvln++)
   {
    h = 0.;   /* will be harmonic mean of kappa in this interval */  

    for (n=0;n<elemn;n++) 
     {
      for (m=m1[n];m<11111;m++) 
        { if (nu[n][m] > nulim[intrvln]) break; }
      m1[n] = m - 1;
      for (m=m1[n];m<11111;m++) 
        { if (nu[n][m] > nulim[intrvln]) break; }
      dm12[n] = m - m1[n] - 1;
     }  /* n-loop */

    for (p=0;p<400;p++)  
     {  
      c = 0.00125 + ((double) p)/400.;  
      nux = (1. - c)*nulim[intrvln] + c*nulim[intrvln+1];    
      g = 0.;   /* will be sum of kappa at this value of _nu_  */

      for (n=0;n<elemn;n++)
       {

        for (m = m1[n];m<=11111;m++) { if (nu[n][m] > nux) break; }
        m1[n] = m - 1;  

        a = (nux - nu[n][m1[n]])/(nu[n][m1[n]+1] - nu[n][m1[n]]);   
        e = (1. - a)*sss[n][m1[n]] + a*sss[n][m1[n]+1];  

        if (dm12[n] > 3)   /* inessential adjustment (?) */
         {
          f = sss[n][m1[n]];
          if (nu[n][m1[n]+1] - nux < nux - nu[n][m1[n]])
             f = sss[n][m1[n]+1];
          e = 0.5*(e + f);  
         }
        g = g + e*dens_atom[n];  

       }   /* n-loop */   
      h = h + 1./g;  

     }  /* p-loop */

    kappa_nu = 400./(h*rho);       /* local harmonic mean of kappa */
    c = (kappa_nu - kappa_sc)/kappa_sc;        /* ratio k_abs/k_sc */    
    nu_intvl = 0.5*(nulim[intrvln] + nulim[intrvln+1]);   

    sprintf(str,"%9.4lf %11.6lf %11.5lf \n",nu_intvl,kappa_nu,c);   
    for (i=0;i<200;i++) if (str[i] == '\0') break;
    write(fdout,str,i);

       /* Now, do the sum for kappa_rm, the Rosseland mean.  
          If  u = h*nu/k*T, the Rosseland mean is given by
          1/kappa_rm = (15/4*pi^4)*{integral of {F(u)/kappa(nu)}*du},
          where F(u) = u^4*exp(u)/{exp(u) - 1}^2.       */
    if (!intrvln)
     { 
      u = ufac*nulim[1]; 
      sum_rm = (u*u*u)/(3.*kappa_nu);    /* crude estimate of IR contrib */ 
     }
    else
     {
      u = ufac*nu_intvl;     /* h*nu/k*T */    
      du = ufac*(nulim[intrvln+1] - nulim[intrvln]);   
      a = exp(u);   b = (a - 1.);  c = u*u; 
      sum_rm = sum_rm + (a*c*c*du)/(b*b*kappa_nu);   
     }

   }  /* intrvln-loop */

  kappa_rm = 25.975758/sum_rm;   /* num factor is 4*pi^4/16  */
  c = (kappa_rm - kappa_sc)/kappa_sc;   /* k_abs/k_sc */ 
  sprintf(str,"kappa_Rosseland:  %.6lf,\nk_abs_rm/k_sc = %.6lf.\n\n---\n",kappa_rm,c); 
  for (i=0;i<200;i++) if (str[i] == '\0') break;
  write(fdout,str,i);


  close(fdout);    
  printf("\n---  Output file is \"%s\"  ---\n",fnameout);  


 }  /* --- end program --- */

/* ----------------------------------------------------------------- */