RBS to GeSi profile conversion with nawk

##############################################################
# RBS_to_SiGe_profile:
# 
#   This program attempts to extract material
# profile information from SiGe structure RBS
# data.
# 
# 31May-2June 1993   Dan Connelly djconnel@leland.stanford.edu
# 17Jan1997: changed to L'Ecuyer scattering cross-section; previous model was
#            flawed.
# last modified: 03Feb1997
# 
# 
# Function notation:
# 
# function fname(x1,x2,....xn,\
#                y1,y2,....ym)
# 
# x1,x2,...xn : expected arguments
# y1,y2,...yn : local variables -- will be ignored if passed
# 
#
##############################################################

BEGIN{
  #step down NR to start it at zero
  NR--;

  # set options
  abort_on_zero_Ge_flag = 0;
}

# skip empty lines
(NF==0){
  NR--;
  next;
}

(NR==0){
  # material array index tags
  Si    = 1;          #tag
  Ge    = 2;          #tag

  # material constants
  pi    = atan2(0,-1);
  M[Si] = 28.08;      #atomic weight of Si in amu -- averaged over isotopes
  M[Ge] = 72.59;      #atomic weight of Ge in amu -- averaged over isotopes
  M_He  =  4;         #mass of alpha particles
  Z[Si] = 14;         #nuclear charge on Si
  Z[Ge] = 32;         #nuclear charge on Ge
  Z_He  =  2;         #nuclear charge on alpha particle beam
  N[Si] = 4.978e10;   #number density, um^-3, from RUMP, 1.660440e-24 g/amu ("units")
  N[Ge] = 4.429e10;
  D[Si] = 2.32;       #calculated from M,N -- g/cm^3
  D[Ge] = 5.338;

  e2 = 1/137;          #square of electron charge, in h-bar c = 0

  # correction factors to electronic stopping power and scattering cross-section
  if (dEdr_Si)
    dEdr_factor[Si] = dEdr_Si;
  else
    dEdr_factor[Si] = 1;
  if (dEdr_Ge)
    dEdr_factor[Ge] = dEdr_Ge;
  else
    dEdr_factor[Ge] = 1;
  if (dsdo_Si)
    dsdo_factor[Si] = dsdo_Si;
  else
    dsdo_factor[Si] = 1;
  if (dsdo_Ge)
    dsdo_factor[Ge] = dsdo_Ge;
  else
    dsdo_factor[Ge] = 1;

  # Temp: The following lines are for temporary use only -- they
  # should be user-settable.
  E_He        = 2.20; #energy of alpha particle beam


  # for tilt spectrum -- detecter, wafer normal, and incoming beam are coplanar
  # constraint: theta_deg = 180 - (theta_out - theta_in)
  if (tiltflag) {
    theta_gauge_deg = tiltangle;   # reading on the guage
    theta_in_deg  = theta;         # incoming beam, off-normal
    theta_out_deg = 181.5 - theta_gauge_deg  + theta_in_deg;  # outgoing beam, off-normal
    theta_deg     = 180 + theta_in_deg - theta_out_deg;   # scattering angle
    }

  # for backscattering data -- detecter, wafer normal, and incoming beam form
  # constraint on angles: cos[theta_out] = cos[theta]_in * cos[theta]
  # right angle
  if (bsflag){
    theta_in_deg  = -6;      # incoming beam, off-normal
    theta_deg     = 170.0;   # outgoing beam, off-normal
    theta_out_deg = 10.7;      # scattering angle
    }

  theta_in      = theta_in_deg  * pi / 180;
  theta_out     = theta_out_deg * pi / 180;
  theta         = theta_deg * pi / 180;
  sec_in        = 1/cos(theta_in);
  sec_out       = 1/cos(theta_out);

 
  if (abs(sin(theta/2))<1e-10) {
    system("echo_stderr \"ERROR -- Backscattering at 180deg not allowed -- aborting.\"");
    exit;
  }

  # calculate theta-independent screening factor
  calculate_dsdo_parameters(Si);
  calculate_dsdo_parameters(Ge);

  # determine column labels
  col_chan = 0;
  col_y_Si = 0;
  col_s_Si = 0;
  col_y_Ge = 0;
  col_s_Ge = 0;

  for (i=0;i++<NF;)
    if ($i ~ "^chan.*")
      col_chan = i;
    else if ($i == "y_Si")
      col_y_Si = i;
    else if ($i == "s_Si")
      col_s_Si = i;
    else if ($i == "y_Ge")
      col_y_Ge = i;
    else if ($i == "s_Ge")
      col_s_Ge = i;

  errorcount = 0;
  if (col_chan == 0) {
    system("echo_stderr \"ERROR : Column for channel number not found.\"");
    errorcount++;
  }
  if (col_y_Si == 0) {
    system("echo_stderr \"ERROR : Column for Si counts not found.\"");
    errorcount++;
  }
  if (col_y_Ge == 0) {
    system("echo_stderr \"ERROR : Column for Ge counts not found.\"");
    errorcount++;
  }
  if (errorcount) {
    system("echo_stderr " sprintf("Fatal error count = %d\n",errorcount));
    system("echo_stderr \"Aborting...\"");
    exit;
  }
  next;
}

#read in data lines
{
  y_Si[$col_chan] = $col_y_Si;
  if (col_s_Si)
    s_Si[$col_chan] = $col_s_Si;
  else
    s_Si[$col_chan] = sqrt(abs($col_y_Si)+1);
  y_Ge[$col_chan] = $col_y_Ge;
  if (col_s_Ge)
    s_Ge[$col_chan] = $col_s_Ge;
  else
    s_Ge[$col_chan] = sqrt(abs($col_y_Ge)+1);
}


# process the data
END{
  # calculate the Si and Ge edge, or used pre-specified values

  cSi = find_surface_channel(y_Si);
  if (chan_si != ""){
    system("echo_stderr \"Si surface channel extraction overridden:\"" chan_si);
    chanmax[Si]=chan_si;
  }
  else{
    if (cSi == 0)
      system("echo_stderr ERROR : Si surface channel not found.");
    chanmax[Si] = cSi;
  }

  cGe = find_surface_channel(y_Ge);
  if (chan_ge != ""){
    system("echo_stderr \"Ge surface channel extraction overridden:\"" chan_ge);
    chanmax[Ge]=chan_ge;
  }
  else{
    if (cGe == 0)
      system("echo_stderr ERROR : Ge surface channel not found.");
    chanmax[Ge] = cGe;
  }

  if ((chanmax[Si] == 0) || (chanmax[Ge] == 0))
    exit;


  # find the starting values
  do_surface_calculations();

  # output the column labels for the output table
  print_column_labels();

  # do the calculations vs depth
  while ((chan[Si]>=0)&&(chan[Ge]>=0)){
    x = add_depth_point();
    if ((abort_on_zero_Ge_flag && (x == 0))||((1*z_max > 0)&&((z_Si>z_max)||(z_Ge>z_max))))
      break;
  }
}

#########################################
# FIND_SURFACE_CHANNEL
#
# Try to determine the channels marking the surface counts
# of Ge, Si.
#

function find_surface_channel(y,\
  xx,yy,dyy,yymax,i,imin,imax,c_edge,dx0,dx1,dx2)
{
  c_edge = 0;
  yymax = "";
  for (i in y) {
    norm = 3;
    yy[i] = 3*y[i];
    xx[i] = 3*i;
    if (y[i-2] != ""){
      yy[i] += y[i-2];
      xx[i] += i-2;
      norm += 1;
    }
    if (y[i-1] != ""){
      yy[i] += 2*y[i-1];
      xx[i] += 2*(i-1);
      norm += 2;
    }
    if (y[i+1] != ""){
      yy[i] += 2*(y[i+1]);
      xx[i] += 2*(i+1);
      norm += 2;
    }
    if (y[i+2] != ""){
      yy[i] += y[i+2];
      xx[i] += i+2;
      norm += 1;
    }
    xx[i]/=norm;
    yy[i]/=norm;

    if ((yymax == "") || (yy[i] > yymax))
      yymax = yy[i];

  }

  # calculate approximation to dy/dx
  for (i in y)
    if ((1*yy[i] != 0) &&\
	(xx[i+1] != "") &&\
	(xx[i-1] != "") &&\
	(xx[i] != "") &&\
	((dx2=xx[i+1]-xx[i])!=0) &&\
	((dx0=xx[i]-xx[i-1])!=0) &&\
	((dx1=xx[i+1]-xx[i-1])!=0))
      dyy[i] = ((yy[i+1]-yy[i])*dx0/dx2 +\
		(yy[i]-yy[i-1])*dx2/dx0)/\
	          dx1;

  imax = "";
  imin = "";
  for (i in y)
    if (y[i] != "") {
      if ((imax == "") || (int(i) > int(imax)))
	imax = int(i);
      if ((imin == "") || (int(i) < int(imin)))
	imin = int(i);
    }

  if (imax > imin+3)
  for (i = imax-3; i >= imin; i--)
    if (dyy[i] != "")
      if (yy[i] > yymax/10)
	if ((dyy[i] > dyy[i+1]) && (dyy[i] > dyy[i+2]) && (dyy[i] > dyy[i+3])) {
	  system("echo_stderr Edge found at channel " (c_edge = int(i+3.5)));
	  break;
	}
    
  return c_edge;
}
       


##########################################
# DO_SURFACE_CALCULATIONS
#
# calculate parameters at surface
#
# Global arrays initialized here:
#   z_array[Ge,0]  : the depth at the surface = 0
#   z_array[Si,0]  : the depth at the surface = 0
#   x_array[Ge,0]  : the alloy content at the surface
#   x_array[Si,0]  : the alloy content at the surface
#   dEindz[Ge,0]   : the electronic stopping rate at the surface
#   dEindz[Si,0]   : the electronic stopping rate at the surface
#
# Global variables initialized here:
#   Echan0    : energy of channel zero
#   Echan1    : energy of channel 1 - Echan0
#   chan[Si]  : channel number of next point under consideration
#               in silicon data
#   chan[Ge]  : channel number of next point under consideration
#               in germanium data
#

function do_surface_calculations(\
				 xmean,xmax,xmin){

  # print comment line to be stripped with stripcomments, for example
  printf("/************************************\n");

  # print simulation parameters
  printf "%-31s= %s amu\n","M[Si]",M[Si];
  printf "%-31s= %s amu\n","M[Ge]",M[Ge];
  printf "%-31s= %s amu\n","M[He]",M_He;
  printf "%-31s= %s\n","Z[Si]",Z[Si];
  printf "%-31s= %s\n","Z[Ge]",Z[Ge];
  printf "%-31s= %s\n","Z[He]",Z_He;
  printf "%-31s= %s g/cm^3\n","Density[Si]",D[Si];
  printf "%-31s= %s g/cm^3\n","Density[Ge]",D[Ge];
  printf "%-31s= %s MeV\n","He energy",E_He;
  printf "%-31s= %s\n","Si surface channel",chanmax[Si];
  printf "%-31s= %s\n","Ge surface channel",chanmax[Ge];
  printf "%-31s= ","Tilt/Backscatter mode";
  if (tiltflag)
    print "tilt";
  else if (bsflag)
    print "backscatter";
  else
    print "???";
  if (tiltflag)
    printf "%-50s= %s deg\n","Angle reading on gauge",theta_gauge_deg;
  printf "%-50s= %s deg\n","Angle of incoming beam to wafer normal",theta_in_deg;
  printf "%-50s= %s deg\n","Angle of detected beam to wafer normal",theta_out_deg;
  printf "%-50s= %s deg\n","Scattering angle: incoming to outgoing beams",theta_deg;

  # print header comments
  printf "%-50s= %s MeV\n","Recoil energy of He off surface Si",\
    Emax[Si] = E_He * (elastic_factor[Si] = elasticity(M[Si],theta));
  printf "%-50s= %s MeV\n","Recoil energy of He off surface Ge",\
    Emax[Ge] = E_He * (elastic_factor[Ge] = elasticity(M[Ge],theta));
  printf "%-31s= %g MeV\n","Difference",Emax[Ge]-Emax[Si];
  printf "%-31s= %s\n","Elasticity He off Si",elastic_factor[Si];
  printf "%-31s= %s\n","Elasticity He off Ge",elastic_factor[Ge];
  printf "%-31s= %g MeV\n","Energy of channel 0",\
    Echan0=(chanmax[Ge]*Emax[Si] - chanmax[Si]*Emax[Ge])/(chanmax[Ge]-chanmax[Si]);
  printf "%-31s= %g MeV\n","E[channel 1] - E[channel0]",
    Echan1=(Emax[Si]-Emax[Ge])/(chanmax[Si]-chanmax[Ge]);

  printf "%-31s= %g\n","dE/dr correction factor for Si",dEdr_factor[Si];
  printf "%-31s= %g\n","dE/dr correction factor for Ge",dEdr_factor[Ge];

  printf("dE/dr in Si at E_He            = %g MeV/um\n", temp0=dEdr(0,E_He));
  printf("dE/dr in Ge at E_He            = %g MeV/um\n", temp1=dEdr(1,E_He));
  printf("ratio dE/dr at E_He Ge/Si      = %g\n", temp1/temp0);
  printf("dE/dr in Si at Si edge         = %g MeV/um\n", temp0=dEdr(0,Emax[Si]));
  printf("dE/dr in Ge at Si edge         = %g MeV/um\n", temp1=dEdr(1,Emax[Si]));
  printf("ratio dE/dr at Si edge Ge/Si   = %g\n", temp1/temp0);
  printf("dE/dr in Si at Ge edge         = %g MeV/um\n", temp0=dEdr(0,Emax[Ge]));
  printf("dE/dr in Ge at Ge edge         = %g MeV/um\n", temp1=dEdr(1,Emax[Ge]));
  printf("ratio dE/dr at Ge edge Ge/Si   = %g\n", temp1/temp0);
  printf("dE/dr in Si at chan 0          = %g MeV/um\n", temp0=dEdr(0,Echan0));
  printf("dE/dr in Ge at chan 0          = %g MeV/um\n", temp1=dEdr(1,Echan0));
  printf("ratio dE/dr at chan 0 Ge/Si    = %g\n", temp1/temp0);

  printf("d Sigma_Si/d Omega at E_He     = %g MeV^-2\n",
	 dsdo(Si,E_He));
  printf("d Sigma_Si/d Omega at chan 0   = %g MeV^-2\n",
	 dsdo(Si,Echan0));
  printf("d Sigma_Ge/d Omega at E_He     = %g MeV^-2\n",
	 dsdo(Ge,E_He));
  printf("d Sigma_Ge/d Omega at chan 0   = %g MeV^-2\n",
	 dsdo(Ge,Echan0));
  printf("d Sigma_Ge/d Sigma_Si at E_He  = %g\n",
	 dsdo(Ge,E_He)/dsdo(Si,E_He));
  printf("d Sigma_Ge/d Sigma_Si at chan 0= %g\n",
	 dsdo(Ge,Echan0)/dsdo(Si,Echan0));


  printf("Count at surface : Si(channel %d)= %d, Ge(channel %d)= %d\n",
	 chanmax[Si],
	 y_Si[chanmax[Si]],
	 chanmax[Ge],
	 y_Ge[chanmax[Ge]]);

  printf("Count 5 channels lower : Si(channel %d)= %d, Ge(channel %d)= %d\n",
	 chanmax[Si]-5,
	 y_Si[chanmax[Si]-5],
	 chanmax[Ge]-5,
	 y_Ge[chanmax[Ge]-5]);

  xmean= find_surface_x(y_Si[chanmax[Si]],y_Ge[chanmax[Ge]]);
  d_array[Si,0] =\
  d_array[Ge,0] = net_density;   /* set by routine */

  xmax = find_surface_x(y_Si[chanmax[Si]]-s_Si[chanmax[Si]],y_Ge[chanmax[Ge]]+s_Ge[chanmax[Ge]]);
  xmin = find_surface_x(y_Si[chanmax[Si]]+s_Si[chanmax[Si]],y_Ge[chanmax[Ge]]-s_Ge[chanmax[Ge]]);

  printf("x at surface       = %g (sqrt2-sigma range= %g,%g)\n",xmean,xmin,xmax);
  printf("sigma of surface x = %g\n",(xmax-xmin)/sqrt(8));

  # print end-of-comment line
  printf("************************************/\n");

  # initialize global arrays
  z_array[Si,0]  = \
  z_array[Ge,0]  = 0;
  x_array[Si,0]  = \
  x_array[Ge,0]  = xmean;
  Ein[Si,0]      = \
  Ein[Ge,0]      = E_He;
  dEindz[Si,0]   = \
  dEindz[Ge,0]   = dEdr(xmean,E_He)*sec_in;

  #initialize global variables
  npoints[Ge]  = 0;        # there are 0 new points in the depth and energy arrays
  npoints[Si]  = 0;        # there are 0 new points in the depth and energy arrays
  chan[Si] = chanmax[Si]   # we looked last at this point in Si data
  chan[Ge] = chanmax[Ge]   # we looked last at this point in Ge data
  ns_last  = 0;            # next Si point to be printed out

  #
  # convert net density to something proportional to atomic density
  #  atoms/volume =
  #        (dNcounts/dEout) * (dEout/dz)                                 1
  #   --------------------------------------------  -----------------------------------------
  #    (dr/dz = sec(theta_in)) * (d sigma/d omega)   (omega_detector) * (N_He in per channel)
  #
  #  the first is the thing we calculate.
  #
  # units: (counts/MeV) (MeV/um)
  #        --------------------- =  um^-3 Str-counts
  #              um^2/Str            
  #
  # Note: for units to work e^2 must be set appropriately and cross-section formula
  #      needs to be adjusted gaussian -> SI

  density_factor = 1/(sec_in*Echan1);

}

########################################
# FIND_SURFACE_X
# this routine tries to determine the value
# of x at the surface
#

function find_surface_x(ySi,yGe,\
			xGe,xSi,zGe,zSi,dedz,x,i,imax) {

  # variables:
  # ----------
  #   ySi,yGe : counts in channel
  #   zSi,zGe : counts in channel over differential scattering cross section
  #   xSi,xGe : zSi,zGe divided by differential change in round-trip energy
  #             change with depth
  #   dedz    : rate of change of energy of beam coming out with respect to
  #             z-axis depth -- note this must be scaled by the |sec| of the
  #             angle to convert from de/d(ion path length)
  # return value: ge mole fraction
  # net_density : number proportional to the observed net particle density --
  #             global variable set by this routine

  # guess at x_surface and then iterate
  zGe = yGe / \
    dsdo(Ge,E_He);
  x = zGe / \
    ((zSi = ySi / \
      dsdo(Si,E_He)) + \
     zGe);

  for (i=(imax=3);i--;) {
    dedz = dEdr(x,E_He)*sec_out;

    # net dEdz where x is the depth at which the backscatter occurs:
    #   (a) electron damping going in the extra dx +
    #   (b) change in inelastic component of scatter at x+dx vs x +
    #   (c) electron damping going out dx from the scatter site +
    #   (d) change in electron going out from x to zero
    # At the surface (d) is zero
    #
    # note the change in the scattering energy, since the elasticity
    # is energy-independent, is (elasticity-1)(dE/dx going in) (less energy
    # lost in the collision the deeper one goes).
    # This adds to dE/dx going in and out to yield
    #
    # net dE/dx = ((elasticity)(dE/dx going in)) + dE/dx going out

    xGe = zGe*(elastic_factor[Ge]*dedz + dEdr(x,Emax[Ge])*sec_in);
    xSi = zSi*(elastic_factor[Si]*dedz + dEdr(x,Emax[Si])*sec_in);
    x = xGe/ \
      net_density=(xSi+xGe);
  }
  return x;
}

#################################################
# ADD_DEPTH_POINT
# This routine extends the array of (x,E) values
# for alphas bouncing off Ge atoms.  It uses the
# array of x versus depth z...
#
# global arrays:
#   z_array[Ge,n]  : depth of nth data point of Ge
#   z_array[Si,n]  : depth of nth data point of Si
#   x_array[Ge,n]  : x value at nth data point of Ge
#   x_array[Si,n]  : x value at nth data point of Si
#   Ein[Ge,n]      : energy of alphas at the grid point for Si 
#   Ein[Si,n]      : energy of alphas at the grid point for Ge
#   dEindz[Ge,n]   : damping rate of incoming alphas at the grid point of Ge
#   dEindz[Si,n]   : damping rate of incoming alphas at the grid point of Si
#
# global variables:
#   npoints[Ge]    : number of points in the depth array and x_array of Ge
#   npoints[Si]    : number of points in the depth array and x_array of Si
#   chan[Ge]       : this is stepped down as a new data point is added
#   chan[Si]       : channel of Si channel under consideration -- channel
#                    is chosen to match germanium channel as closely as
#                    possible
#
# local arrays
#   Eout[Si or Ge,n]   : energy of alphas at the grid point on their way out
#                        assuming they reach the surface at the energy of interest
#   dEoutdz[Si or Ge,n]: rate of energy loss of alphas on their way out assuming
#                        they reach the surface at the energy of interest
#

function add_depth_point(\
                         i,n,delta_z,x,E,Eout,dEoutdz,dedz,ns,ng,deltaz1,deltaz2,deltaz,zg1,zg2,zs1,zs2) {

  # step down Ge channel number and determine the energy of that channel
  E = channel_to_energy( \
      --chan[Ge] \
      );

  if (chan[Ge]<0)
    return 0;

  Eout[Ge,0]    = E;                             # particle reaches surface at chosen energy from Ge impact
  dEoutdz[Ge,0] = dEdr(x_array[Ge,0],E)*sec_out; # rate at which particle is losing energy at the surface

  for(n=0;n++<npoints[Ge];) {
    delta_z   = z_array[Ge,n]-z_array[Ge,n-1];          # change in depth over the interval

    dE = delta_z*dEoutdz[Ge,n-1];
    for(i=3;i--;){
      dedz = dEdr(x_array[Ge,n],Eout[Ge,n-1]+dE)*sec_out;
      dE = delta_z * (dedz + dEoutdz[Ge,n-1])/2;
    }
    dEoutdz[Ge,n] = dedz;
    Eout[Ge,n] = Eout[Ge,n-1] + dE;
  }
  n--;

  # extrapolate x value of new point --
  # we want Eout(z) = Ein(z) * elasticity
  #
  # approximate
  #   Eout  = elasticity Ein = Eout0 + delta_z dEout/dz
  #                      Ein = Ein0  - delta_z dEin/dz
  #
  # substitute for Ein: 
  # elasticity (Ein0 - delta_z dEin/dz) = Eout0 + delta_z dEout/dz
  #
  # combine terms:
  # delta_z (elasticity dEin/dz + dEout/dz) = elasticity Ein0 - Eout0
  #
  # divide:
  # delta_z = (elasticity Ein0 - Eout0) / (elasticity dEin/dz + dEout/dz) 

  z_array[Ge,n+1] = z_array[Ge,n] + \
    (elastic_factor[Ge]*Ein[Ge,n] - Eout[Ge,n]) / (elastic_factor[Ge]*dEindz[Ge,n] + dEoutdz[Ge,n]);

  # step until we get to a silicon point past the last germanum one.
  while(z_array[Si,npoints[Si]] <= z_array[Ge,npoints[Ge]]) {

    # if we run out of silicon data call it quits, otherwise step down silicon grid
    E = channel_to_energy( \
        --chan[Si]\
        );

    if (chan[Si] < 0)
      return 0;

    Eout[Si,0]    = E;                             # particle reaches surface at chosen energy from Ge impact
    dEoutdz[Si,0] = dEdr(x_array[Si,0],E)*sec_out; # rate at which particle is losing energy at the surface

    for(n=0;n++<npoints[Si];) {
      delta_z = z_array[Si,n]-z_array[Si,n-1];        # change in depth over the interval
	  
      dE = delta_z*dEoutdz[Si,n-1];

      for(i=3;i--;){
	dedz = dEdr(x_array[Si,n],Eout[Si,n-1]+dE)*sec_out;
	dE = delta_z * (dedz + dEoutdz[Si,n-1])/2;
      }
      dEoutdz[Si,n] = dedz;
      Eout[Si,n] = Eout[Si,n-1] + dE;
    }
    n--;
    ns = ++npoints[Si];
    ng = npoints[Ge];

    # extrapolate the depth, as with the germanium data
    zs2 = \
      z_array[Si,ns] = \
	(zs1 = z_array[Si,n]) + \
	  (elastic_factor[Si]*Ein[Si,n] - Eout[Si,n]) / (elastic_factor[Si]*dEindz[Si,n] + dEoutdz[Si,n]);

    if (ng){
      zg1 = z_array[Ge,ng-1];
      zg2 = z_array[Ge,ng];

      # if it exists, calc d_array for straddled Si point
      if (ns>1){
        # interpolate calc net density for straddled Si point
        deltaz = \
	  (deltaz1 = (zs1 - zg1)) + \
	    (deltaz2 = (zg2 - zs1));

	  d_array[Si,n] = \
	    (deltaz2 * d_array[Ge,ng-1] + \
	     deltaz1 * d_array[Ge,ng]) / \
	       deltaz;
      }

      # determine x and Ein for the new Si point

      deltaz = \
	(deltaz1 = (zs2 - zg1)) + \
	  (deltaz2 = (zg2 - zs2));
      
      x_array[Si,ns] = \
	(deltaz2 * x_array[Ge,ng-1] + \
	 deltaz1 * x_array[Ge,ng]) / \
	   deltaz;

      # don't bother calculating d_array[Si,ns], which isn't used.

      Ein[Si,ns]=Ein[Si,n]+dEindz[Si,n]*(zs2-zs1);
      for (i=2;i--;){
	dEindz[Si,ns] = sec_in*dEdr(x_array[Si,ns],Ein[Si,ns]);
	Ein[Si,ns] = Ein[Si,n] - (dEindz[Si,n]+dEindz[Si,ns])*(zs2-zs1)/2;
      }
    }
    else {
      dEindz[Si,ns] = dEindz[Si,ns-1];
      for (i=2;i--;)
	dEindz[Si,ns] = sec_in * dEdr(\
	  (x_array[Si,ns] = x_array[Si,0]),\
          (Ein[Si,ns] = Ein[Si,ns-1] +(z_array[Si,ns-1]-z_array[Si,ns])*(dEindz[Si,ns-1]+dEindz[Si,ns])/2)\
	  );
    }

  }

  #################################  
  # Now we have the last two Si points straddling the last Ge point.  We can extrapolate
  # the Si counts to match the depth and calculate the atomic germanium fraction.
  # we can do the same for the second-to-last Si point using the straddling germanium
  # data
  #
  # variables:
  # ----------
  #   ySi,yGe : counts in channel
  #   zSi,zGe : counts in channel over differential scattering cross section
  #   xSi,xGe : zSi,zGe divided by differential change in round-trip energy
  #             change with depth
  

  #############################################################
  # Calculate results for new previous then for current Ge point

  ns = npoints[Si];

  for (ng=max(npoints[Ge],1);ng<=npoints[Ge]+1;ng++) {
    dE = \
      (dEindz[Ge,ng] = dEindz[Ge,ng-1]) * \
	(dz_Ge = z_array[Ge,ng]-z_array[Ge,ng-1]);
    
    for (i=2;i--;){
      #calculate Ein
      Ein[Ge,ng]= Ein[Ge,ng-1] - dE;
      
      deltaz = \
	(deltaz1 = (z_array[Ge,ng] - z_array[Si,ns-1])) +\
	  (deltaz2 = (z_array[Si,ns] - z_array[Ge,ng]));
      
      zGe = \
	(yGe = y_Ge[chanmax[Ge]-ng]) /\
	  dsdo(Ge,Ein[Ge,ng]);
      zSi = \
	(ySi = ((y_Si[chanmax[Si]-ns]  * deltaz1 + \
		 y_Si[chanmax[Si]-ns+1]* deltaz2)/ \
		deltaz)) / \
		  dsdo(Si, \
		       (EinSi = \
			(Ein[Si,ns]   * deltaz1 + \
			 Ein[Si,ns-1] * deltaz2)/ \
			deltaz));
      
      # numerically determine dedz
      dedz_Si = Echan1 / (z_array[Si,ns] - z_array[Si,ns-1]);
      dedz_Ge = Echan1 / (z_array[Ge,ng] - z_array[Ge,ng-1]);

      x = x_array[Ge,ng] = zGe*dedz_Ge / \
	d_array[Ge,ng] = (zGe*dedz_Ge + zSi*dedz_Si);
      dEindz[Ge,ng] = dEdr(x,Ein[Ge,ng])*sec_in;
      dE = (dEindz[Ge,ng-1] + dEindz[Ge,ng]) * dz_Ge / 2;
    }

    # revise all Si points affected by the new Ge point
    for (n=ns; n && (z_array[Si,n]>z_array[Ge,ng-1]); n--);
    while(n++<ns) {
      zs  = z_array[Si,n];
      zg1 = z_array[Ge,ng-1];
      zg2 = z_array[Ge,ng];
      
      # interpolate x and calc Ein for the straddled Si point based on new Ge point data
      deltaz = \
	(deltaz1 = (zs - zg1)) + \
	  (deltaz2 = (zg2 - zs));
      
      x_array[Si,n] = \
	(deltaz2 * x_array[Ge,ng-1] + \
	 deltaz1 * x_array[Ge,ng]) / \
	   deltaz;
	   
      for (i=2;i--;){
	dEindz[Si,n] = sec_in*dEdr(x_array[Si,n],Ein[Si,n]);
	Ein[Si,n] = Ein[Si,n-1] - (dEindz[Si,n-1]+dEindz[Si,n])*(z_array[Si,n]-z_array[Si,n-1])/2;
      }
    }
    
    

  }
  npoints[Ge]++;
  ng--;


  #######################################
  # output the data point
  #
  i=chanmax[Ge]-ng+1;
  print "Ge",\
    i,\
    ng-1,
    y_Si[i],\
    y_Ge[i],\
    y_Ge[i]+y_Si[i],\
    (z_Ge = z_array[Ge,ng-1]),\
    (x_last = x_array[Ge,ng-1]),\
    Ein[Ge,ng-1],\
    Eout[Ge,0],\
    dEindz[Ge,ng-1],\
    dedz_Ge,\
    d_array[Ge,ng-1]*density_factor;

  for (n=ns_last; n<ns; n++) {
    i=chanmax[Si]-ns+1;
    print "Si",\
      i,\
      ns-1,
      y_Si[i],\
      y_Ge[i],\
      y_Ge[i]+y_Si[i],\
      (z_Si = z_array[Si,n]),\
      (x_last = x_array[Si,n]),\
      Ein[Si,n],\
      Eout[Si,0],\
      dEindz[Si,n],\
      dedz_Si,\
      d_array[Si,n]*density_factor;

   }
  ns_last = ns;
 
  return x_last;

}

##################################################
# PRINT_COLUMN_LABELS
#
function print_column_labels(){
print\
  "atom",\
  "chan",\
  "chan_rel",\
  "counts_Si",\
  "counts_Ge",\
  "counts",\
  "depth",\
  "x",\
  "Escatter",\
  "Eout",\
  "dEscatter/dz",\
  "dEout/dz",\
  "net_density"
}

##################################################
# SQUARE : x*x
#
function square(x){
  return x*x;
}

##################################################
# ABS
#
function abs(x){
  if (x<0)
    return -x;
  else
    return x;
}

##################################################
# MIN
#
function max(x,y){
  if (x<y)
    return y;
  else
    return x;
}

##################################################
# CHANNEL_TO_ENERGY : convert channel number
# to energy.  Expect there is a linear relation
# and the coefficients have been set
#
function channel_to_energy(channel){
  return Echan0 + channel*Echan1;
}

##################################################
# DSDO : calculate the differential cross section of
# the particular element to alpha particles at the given
# angle for the given target mass at the given energy
# with the given nuclear charge of the target.
#
function dsdo(species,E,\
              xcom,xlab,Ecom){
  xcom = (e2 * Z_He * Z[species] / (4 * E * stheta2))^2;
  xlab = xcom * (r1mr2s2 + ctheta)^2/r1mr2s2;
  Ecom = E*Efactor
  return (xlab * (1 - ((Z[species]*Z_He)^(4/3))*49e-6/Ecom));
  
}
#################################################
# CALCULATE_DSDO_PARAMETERS:
# set global variables appropriate for calculation
# of scattering cross-section
#
# These parameters are based on the L'Ecuyer
# lab-frame approximation to the screening-corrected
# Rutherford Backscattering formula
#
function calculate_dsdo_parameters(species){
  rho     = M_He/M[species];
  Efactor = (1 + rho)^-2;
  stheta2 = sin(theta)^2;
  ctheta  = cos(theta);
  r1mr2s2 = sqrt(1 - rho^2*stheta2);
}

#################################################
# ELASTICITY
# return the collision elasticity factor for a given
# target versus alpha particles\
#
function elasticity(M,theta,\
                    Mcos,Msum){
  Mcos = M_He*cos(theta) / \
   (Msum = M + M_He);

  return square(Mcos + sqrt(Mcos*Mcos + (M-M_He)/Msum));
}
  

##################################################
# STOPPING_POWER : return the stopping power
# at the given energy.
#

function stopping_power(species,e,\
                        slow, shigh){

  if (species==Ge) {
   slow  = 143.67 * e^0.466;
   shigh = (79.24/e)*log(1 + (1.185/e) + (7.993*e));
  }
  else if (species==Si) {
   slow  = 187.16 * e^0.65;
   shigh = (49.34/e)*log(1 + (1.788/e) + (4.133*e));
  }
  else {
   print "ERROR : STOPPING POWER CALLED WITH ILLEGAL SPECIES";
   return 0;
  }

  return 10^-13*N[species]/(1/slow + 1/shigh);
}

##################################################
# dEdr: Find net stopping power for alloy of a given
# composition at the given energy assuming a linear
# superposition of stopping powers.
#
# Note : this returns the change in E with respect to
#    r, the path length, not z, the depth.  Multiply
#    by abs(sec(theta_{in,out})) to get dEdz.
#
function dEdr(x,E) {

  # assume Vegards law applies to number density instead of
  # lattice constant.

  return \
    x     * stopping_power(Ge,E) * dEdr_factor[Ge] + \
    (1-x) * stopping_power(Si,E) * dEdr_factor[Si];
}