function get_rcoefs, ener, graz_angle, sigma, comp, params, cathodemod, $
                     interface=interface, debug=debug, dump=dump, nroughcor=nroughcor
;
;  Makes a call of 'mk_recipe' and the params arguments are
;  described there
;
;     comp is composition, 0 for Pt/C, 1 for W/Si, and 2 for NiV/C
;
;     nroughcor returns 1 if it has been necessary to correct for roughgrad <= rosub
;

   common optical_constants, cimd_c_real, cimd_c_imag, $
                             cimd_pt_real, cimd_pt_imag, $
                             cimd_si_real, cimd_si_imag, $
                             cimd_w_real, cimd_w_imag, $
                             cimd_niv_real, cimd_niv_imag, $
                             cimd_sicrystal_real, cimd_sicrystal_imag, $
                             cimd_escale

   if( n_params() lt 5 ) then begin
      print,'Syntax: coefs = get_rcoefs( ener, graz_angle, sigma, comp, params[ , cathodemod] )
      print,'         keywords: interface, debug, dump, nroughcor'
      print,'  graz_angle must be in degrees'
      print,'  Optional ''cathodemod'' are modification factors for layer thicknesses'
      print,''
      print,'  ''sigma'' can be a scalar which is then used for all interfaces in the multilayer'
      print,'  If ''sigma'' is an array it must have three or four elements as:'
      print,'       [sigma_HZ, sigma_LZ, sigma_SUBSTR]'
      print,'   or  [sigma_HZ, sigma_LZ, sigma_SUBSTR, roughgrad]'
      print,''
      print,'  ''comp'' 0 - Pt/C,   1 - W/Si,   and 2 - NiV/C'
      print,''
      print,'   The default value for interface is 0 (zero) i.e. error function interface shape'
      print,'   If given it should either be a scalar or an array with N+1 elements of integers'
      print,''
      print,'   The optical constant file names are assumed to be placed in environment variables'
      print,'    C_OPTCSTFILE, PT_OPTCSTFILE, SI_OPTCSTFILE, and W_OPTCSTFILE.'
      return, 0
   endif

   dbg = 0
   if keyword_set(debug) then dbg = 1
   dmp = 0
   if keyword_set(dump) then begin
      dmp = 1
      openw,lun,'dump_get_rcoefs.txt',/get_lun
   endif

   if not keyword_set(interface) then interface = 0

   nroughcor = 0 ;  Flag for roughness adaption to avoid NaN results (see later)

   nan_reported = 0;  Flag for telling if the appearence of one or more NaNs already
                   ;  has been reported

;---------------------- Initialize optical_constants, but only the first time

   if n_elements(cimd_escale) eq 0 then begin
      cimd_escale = xlingen(500, 0.5, 100.)
      cimd_wlscale = 12.3985 / cimd_escale
      ; Use imd_nk to read data

      cimd = imd_nk( getenv('C_OPTCSTFILE'), cimd_wlscale );  Returns dcomplex values
      ;;cimd = imd_nk( 'C_NIST_CPJ_071125', cimd_wlscale );  Returns dcomplex values
      cimd_c_real = double(cimd)
      cimd_c_imag = imaginary(cimd)

      cimd = imd_nk( getenv('SI_OPTCSTFILE'), cimd_wlscale );  Returns dcomplex values
      ;;cimd = imd_nk( 'Si_NIST_CPJ_071125', cimd_wlscale );  Returns dcomplex values
      cimd_si_real = double(cimd)
      cimd_si_imag = imaginary(cimd)

      cimd = imd_nk( getenv('PT_OPTCSTFILE'), cimd_wlscale );  Returns dcomplex values
      ;;cimd = imd_nk( 'Pt_NIST_CPJ_071010', cimd_wlscale );  Returns dcomplex values
      cimd_pt_real = double(cimd)
      cimd_pt_imag = imaginary(cimd)

      cimd = imd_nk( getenv('W_OPTCSTFILE'), cimd_wlscale );  Returns dcomplex values
      ;;cimd = imd_nk( 'W_NIST_CPJ_071010', cimd_wlscale );  Returns dcomplex values
      cimd_w_real = double(cimd)
      cimd_w_imag = imaginary(cimd)

      cimd = imd_nk( 'Ni.93V.07_NIST_CPJ_081013_0.1keV', cimd_wlscale );  Returns dcomplex values
      cimd_niv_real = double(cimd)
      cimd_niv_imag = imaginary(cimd)

      cimd = imd_nk( 'Si_crystal_NIST_CPJ_081013_0.1keV', cimd_wlscale );  Returns dcomplex values
      cimd_sicrystal_real = double(cimd)
      cimd_sicrystal_imag = imaginary(cimd)
   endif

;-----    End of filling common 'optical_constants' ---------

   if dbg then print,'*get_rcoefs* going to call mk_recipe'
   z_arr = mk_recipe( params, debug=debug, dump=dump )
   if dbg then print,'*get_rcoefs* returned from mk_recipe'
   n_layers = n_elements(z_arr)
   n_bilayers = n_layers/2

   ;  ---   Modify the layers if cathode dropouts have occurred

   if n_params() ge 6 then begin
      if dbg then print,'*get_rcoefs* cathode modification'
      if n_elements(cathodemod) lt n_layers then begin
         print,'Bad dimensionality of ''cathodemod'''
         return, 0
      endif
      z_arr = z_arr * cathodemod(0:n_layers-1)
   endif

   theta = 90. - graz_angle;  Must be in degrees

   n_ener = n_elements(ener)

   nc_arr = dcomplexarr(n_layers+2,n_ener)

   nc_arr(0,*) = dcomplex(1.,0.) ; setting the vacuum optical constant

;   Read the nk data and check with 'imd_nk'

   wener = 12.3985 / ener

;------   Light Z layers  ---------

   case comp of
     0 : begin ; Using Pt/C
            elname = 'C'
            creal = interpol( cimd_c_real, cimd_escale, ener )
            cimag = interpol( cimd_c_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
     1 : begin ; Using W/Si
            elname = 'Si'
            creal = interpol( cimd_si_real, cimd_escale, ener )
            cimag = interpol( cimd_si_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
     2 : begin ; Using NiV/C
            elname = 'C'
            creal = interpol( cimd_c_real, cimd_escale, ener )
            cimag = interpol( cimd_c_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
     else : begin
              print,'Illegal composition'
              return, 0
            end
   endcase

   for j = 1, n_layers-1, 2 do nc_arr(j,*) = cimd

;------   Heavy Z layers  ---------

   case comp of
     0 : begin ; Using Pt/C
            elname = 'Pt'
            creal = interpol( cimd_pt_real, cimd_escale, ener )
            cimag = interpol( cimd_pt_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
     1 : begin; Using W/Si
            elname = 'W'
            creal = interpol( cimd_w_real, cimd_escale, ener )
            cimag = interpol( cimd_w_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
     2 : begin; Using NiV/C
            elname = 'NiV'
            creal = interpol( cimd_niv_real, cimd_escale, ener )
            cimag = interpol( cimd_niv_imag, cimd_escale, ener )
            cimd = dcomplex(creal,cimag)
         end
   endcase

   for j = 2, n_layers, 2 do nc_arr(j,*) = cimd

;------   Substrate  ---------

   creal = interpol( cimd_sicrystal_real, cimd_escale, ener )
   cimag = interpol( cimd_sicrystal_imag, cimd_escale, ener )
   cimd = dcomplex(creal,cimag)

   nc_arr(n_layers+1,*) = cimd

   hdr = replicate('//',10)

   k = 0
   hdr(k) = hdr(k) + ' ' + ndate(3);
   k = k + 1
   hdr(k) = hdr(k) + string(theta(0),format='(" theta = ",f8.4," deg")')
   k = k + 1
   hdr(k) = hdr(k) + string(n_ener,ener(0),ener(n_ener-1),$
        format='(i4," energies from ",f8.4," thru ",f8.4," keV")')
   k = k + 1

   ;  Setting the roughnesses

   if n_elements(sigma) gt 2 then begin
      rough_arr = replicate(float(1),n_layers+1)
      rough_arr(0:*:2) = sigma(1)
      rough_arr(1:*:2) = sigma(0)
      rough_arr(n_layers) = sigma(2)
      wnan = where( rough_arr ne rough_arr, nwnan )
      wf = where( 1 - finite( rough_arr ), nwf )
      if nwnan gt 0 or nwf gt 0 then begin
         print,'##13## GET_RCOEFS NaN'
         print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf
         snapshot,v1=sigma,t1='Sigma'
      endif
      if n_elements(sigma) gt 3 then begin
         rosub = 0.5*(sigma(0)+sigma(1))
         c = params(3)
         N = params(0)
         roughgrad = sigma(3)

         ;  Precaution for bad match of rosub and roughgrad

         if roughgrad le rosub then begin
            rosub = 0.95*roughgrad
            nroughcor = 1
         endif

         krough=(rosub/roughgrad)^(1./c)
         brough=(1.-N*krough)/(krough-1.)
         arough=rosub*(brough+N)^c

;   --------   The following section is only for debugging

         for i = 1, fix(N+0.5) do begin
            rq = arough/(brough+i)^c
            nwnan = ( rq ne rq )
            nwf = 1 - finite( rq )
            if nwnan or nwf then begin
               print,'##14## GET_RCOEFS NaN'
               print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf
               snapshot,v1=[rosub,roughgrad,c,N],t1='rosub,roughgrad,c,N', $
                   v2=[krough,brough,arough],t2='krough,brough,arough', $
                   v3=i,t3='index i',v4=rq,t4=' ** rq **'
            endif
         endfor

;   ----------    End of debugging  ----------------

         for i = 1, fix(N+0.5) do rough_arr(2*(i-1):2*i-1) = arough/(brough+i)^c
         wnan = where( rough_arr ne rough_arr, nwnan )
         wf = where( 1 - finite( rough_arr ), nwf )
         if nwnan gt 0 or nwf gt 0 then begin
            print,'##15## GET_RCOEFS NaN'
            print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf
            snapshot,v1=[rosub,roughgrad,c,N],t1='rosub,roughgrad,c,N', $
                v2=[krough,brough,arough],t2='krough,brough,arough'
            nan_reported = 1
         endif
      endif
      hdr(k) = hdr(k) + string(sigma(0),format='(" sigmaHZ = ",f8.4," A")')
      k = k + 1
      hdr(k) = hdr(k) + string(sigma(1),format='(" sigmaLZ = ",f8.4," A")')
      k = k + 1
      hdr(k) = hdr(k) + string(sigma(2),format='(" sigmaSUBSTR = ",f8.4," A")')
   endif else begin
      hdr(k) = hdr(k) + string(sigma,format='(" sigma(common) = ",f8.4," A")')
      rough_arr = sigma
   endelse

   if dmp then begin
      for i = 0, 5 do begin
         printf,lun,wener(i),format='("Lambda",f10.5," A")'
         for j = 0, 5 do begin
            printf,lun,float(nc_arr(j,i)),imaginary(nc_arr(j,i)),format='("        ",f12.9,f12.9)'
         endfor
      endfor
      free_lun,lun
   endif

   ;  ****   The crucial call of the FRESNEL function by Dave Windt:

   fresnel, theta, wener, nc_arr, z_arr, rough_arr, interface, mfc_model=1, ra=out_arr

   wnan = where( out_arr ne out_arr, nwnan )
   wf = where( 1 - finite( out_arr ), nwf )
   if nwnan gt 0 or nwf gt 0 then begin
      print,'##16## GET_RCOEFS NaN'
      print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf
   endif

   return, out_arr

end