libcpixe/stop96.c

Go to the documentation of this file.

/***************************************************************************
    Copyright (C) 2007 by Carlos Pascual-Izarra
    carlos.pascual_AT_users.sourceforge.net                                                  *

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

    -----------------------------------------------------------------
    STOP96.c: Ion stopping subroutines used by hotstop
    By Carlos Pascual-Izarra
    This is a C version of the STOP96 subroutine by J.F.Ziegler for
    being used in TRIM / SRIM. The original from Ziegler was written
    in BASIC.

    IMPORTANT:
    It must be clear that J.F.Ziegler or any other person related to
    TRIM/SRIM have absolutely no responsability of any limitation or
    bug related to the present code.

    The following is the license notice of SRIM:

    "SRIM.EXE, (C) 1984-2002, James F. Ziegler
    Permission to use, copy, modify and distribute
    this software for any non-commercial purpose and
    without fee is hereby granted, provided that
    this copyright and permission notice appear
    on all copies of the software. The name of
    the author may not be used in any advertising
    or publicity pertaining to the use of the software.
    The author makes no warranty or representations
    about the suitability of the software for any
    purpose. It is provided "AS IS" without any express
    or implied warranty, including the implied warranties
    of merchantability, fitness for a particular purpose
    and non-infringement. The author shall not be liable
    for any direct, indirect, special or consequential
    damages resulting from the loss of use, data or
    projects, whether in an action of contract or tort,
    arising out of or in connection with the use or
    performance of this software. Downloading any
    software from this site is implicit acceptance of
    these conditions. The term "SRIM" as used to
    refer to topics concerning the  Stopping and
    Range of Ions in Matter is copyright registered
    in the United States Copyright Office and may
    not be used for commercial purposes without permission.
    For commercial applications, licenses may be obtained
    for incorporating SRIM into applications by writing to
    SRIM.org."

 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>     /* REMEMBER: use -lm when compiling */
#include "stop96.h"  /*Contains the definitions for the stop96.c functions*/
#include "compilopt.h"


void readstopcoef(char *afilename, char *bfilename)
{
/*Fills the scoef[94,55] array with the data in SCOEF.95A and SCOEF.95B files*/
/* scoef[94,55] must be previously declared as a global variable */
/*Note: the row 0 and column 0 of scoef are not used to follow the BASIC standard of labelling arrays: SCOEF(1:93,1:55)*/

FILE *afile, *bfile;
int i,j;
char stemp[300];

/*Open files*/
if (!(afile=fopen(afilename, "r"))){fprintf(stderr, "\ncannot open '%s' for read\n",afilename);exit(1);}
if (!(bfile=fopen(bfilename, "r"))){fprintf(stderr, "\ncannot open '%s' for read\n",bfilename);exit(1);}

/*Skip headers*/
for(i=0;i<2;i++) fgets(stemp,300,afile) ;  /*Skip Headers in file A*/
for(i=0;i<2;i++) fgets(stemp,300,bfile) ;  /*Skip Headers in file B*/

/*fill scoef*/
for(i=1;i<=93;i++) for(j=1;j<=16;j++)fscanf(afile,"%f",&scoef[i][j]);
for(i=1;i<=93;i++) for(j=17;j<=54;j++)fscanf(bfile,"%f",&scoef[i][j]);


/*Close files*/
fclose(afile);
fclose(bfile);
}

float pstop(int Z2, float E)
{

/* scoef[94,55] must be previously declared as a global variable */
float X,SL,SH,ppower,Sproton;
float PE0;
PE0=10;
/*High (E>10000keV/amu) energy stopping*/
if(E>1e4){
  X=logf(E)/E;
  Sproton=scoef[Z2][17]+(scoef[Z2][18]*X)+(scoef[Z2][19]*X*X)+(scoef[Z2][20]/X);
  }
/*Medium energies: (PE0<E<1e4)  , PE0=is the low limit for "medium" energy range (in kev/amu)*/
else if (E>PE0){
  SL=scoef[Z2][9]*powf(E,scoef[Z2][10])+scoef[Z2][11]*powf(E,scoef[Z2][12]);
  SH=scoef[Z2][13]/powf(E,scoef[Z2][14])*logf((scoef[Z2][15]/E)+scoef[Z2][16]*E);
  Sproton=SL*SH/(SL+SH);
  }
/*Low energies (E<PE0) : approx velocity proportional stopping power, conveniently scaled for continuity with medium E range*/
else if(E>0){
  if(Z2<7)ppower=0.35; else ppower=0.45;
  SL=scoef[Z2][9]*powf(PE0,scoef[Z2][10])+scoef[Z2][11]*powf(PE0,scoef[Z2][12]);
  SH=scoef[Z2][13]/powf(PE0,scoef[Z2][14])*logf((scoef[Z2][15]/PE0)+scoef[Z2][16]*PE0);
  Sproton=(SL*SH/(SL+SH))*powf(E/PE0,ppower);
  }
else Sproton=0.;

return(Sproton);
}


float pstop_fit(int Z2,float E,float P0,float P1,float P2,float P3,float P4,float P5,float P6,float P7,float HE0,float HE1,float HE2,float HE3)
{
/*This function is identical to pstop but it gets all the parameters from its arguments instead of looking at scoef[Z2][x]*/
float X,SL,SH,ppower,Sproton;
float PE0;
PE0=10;
/*High (E>10000keV/amu) energy stopping*/
if(E>1e4){
  X=logf(E)/E;
  Sproton=HE0+(HE1*X)+(HE2*X*X)+(HE3/X);
  }
/*Medium energies: (PE0<E<1e4)  , PE0=is the low limit for "medium" energy range (in kev/amu)*/
else if (E>PE0){
  SL=P0*powf(E,P1)+P2*powf(E,P3);
  SH=P4/powf(E,P5)*logf((P6/E)+P7*E);
  Sproton=SL*SH/(SL+SH);
  }
/*Low energies (E<PE0) : approx velocity proportional stopping power, conveniently scaled for continuity with medium E range*/
else if(E>0){
  if(Z2<7)ppower=0.35; else ppower=0.45;
  SL=P0*powf(PE0,P1)+P2*powf(PE0,P3);
  SH=P4/powf(PE0,P5)*logf((P6/PE0)+P7*PE0);
  Sproton=(SL*SH/(SL+SH))*powf(E/PE0,ppower);
  }
else Sproton=0.;

return(Sproton);
}

void stop96(int Z1, int Z2,float Ionweight, float *Ener, int nEner, float *SE){

/*
Z1=incident ion atomic number
Z2=target atomic number
M1=Weight of incident ion (if M1<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
nEner=dim of Ener as well as of SE
SE=pointer to a vector of nEner components where the resultin stoppings will be stored (ev/1e15atcm-2)


*/
float M1,tempf,VFERMI,E,HE0,HE,B,HEH,YRMIN,VRMIN,V,VR,YR,Q,LAMBDA0,LAMBDA1,L,ZETA0,ZETA,VMIN,EEE,Sproton,EION,VFCORR0,VFCORR1,HIPOWER;
int i,j;

/*Do checks*/
if(Z1<1 || Z1>92) {fprintf(stderr, "\nError: Z1 must be between 1 and 92 \n");exit(1);}
if(Z2<1 || Z2>92) {fprintf(stderr, "\nError: Z2 must be between 1 and 92 \n");exit(1);}
if(nEner<1){fprintf(stderr, "\nError: I need at least one energy  \n");exit(1);}


/*Fill some variables with data*/
VFERMI= scoef[Z2][7];       /*VFermi of Target*/
if(Ionweight<1)M1=scoef[Z1][3]; else M1=Ionweight;    /*MOI weight if M1 wasn't specified*/

/*Loop for stopping values*/
for(i=0;i<nEner;i++){
  E=Ener[i]/M1;

  if(Z1==1){ /*stop_H*/
    SE[i]=pstop(Z2,E);
    }


  else if(Z1==2){ /*stop_He */
    HE0=1;                            /*Minimum limit for "medium" energy range for He (units=keV/amu)*/
    if(HE0>E)HE=HE0;else HE=E;        /*HE=max{E,HE0}*/
    B=logf(HE);
    tempf=.2865+.1266*B-.001429*B*B+.02402*B*B*B-.01135*powf(B,4)+.001475*powf(B,5);
    if(tempf>30) tempf=30;
    HEH=1.-expf(-tempf);
    /* ADD Z1^3 EFFECT TO HE/H STOPPING POWER RATIO, HEH */
    if(HE<1)HE=1;  /*I think this assertion is innecessary!!!! (it its in the original code but seems redundant)???????*/
    tempf=(1.+(.007+.00005*Z2)*expf(-powf(7.6-logf(HE),2)));
    HEH*=(tempf*tempf);
    SE[i]=pstop(Z2,HE)*HEH*4.;
    /* If E was lower than HE0 , CALC. HE VELOCITY PROPORTIONAL STOPPING (low energy range)*/
    if(HE0>E)SE[i]*=sqrtf(E/HE0);
    }


  else if(Z1>2){ /*stop_HI*/
    /*USE VELOCITY STOPPING FOR (YRMIN=VR/Z1**.67) <= 0.13  _OR_  VR<=1.0*/
    YRMIN=0.13;   /*Minimum limit for YR*/
    VRMIN=1.0;    /*Minimum limit for VR*/
    V=sqrtf(E/25.)/VFERMI;  /*V is the "relative velocity" but...relative to what????? (Fermi, Bohr,...????)*/
    if(V<1) VR=(3.*VFERMI/4.)*(1+(2.*V*V/3.)-powf(V,4.)/15.);
    else VR=V*VFERMI*(1.+1./(5.*V*V));
    /*put a minimum limit to VR*/
    if(VR<VRMIN)VR=VRMIN;
    /*calc YR as VR/Z1^(2/3) and put a minimum limit to it*/
    YR=VR/powf(Z1,.6667);
    if(YR<YRMIN)YR=YRMIN;
    /*Calculate the Ionization level (Q) of the ion at velocity YR*/
    tempf=-.803*powf(YR,0.3)+1.3167*powf(YR,0.6)+.38157*YR+.008983*YR*YR;
    if(tempf>50.)Q=1; else if(tempf<=0)Q=0; else  Q=1.-expf(-tempf);    /*It prevents underflow and wasting of time  in the Q calc. */
    /*NOW WE CONVERT IONIZATION LEVEL TO EFFECTIVE CHARGE
    ---------------------- Screening Distance of Ion (Lambda in B.& K.)*/
    for(j=22; j<=39 && scoef[93][j]<Q ;j++);j-=1;  /*Find Q Interpolation in SCOEF(93,22-39)*/
    if(j<22)j=22;else if(j>38)j=38; /*Boilerplate ??*/
    LAMBDA0=scoef[Z1][j];
    LAMBDA1=(Q-scoef[93][j])*(scoef[Z1][j+1]-scoef[Z1][j])/(scoef[93][j+1]-scoef[93][j]);
    L=(LAMBDA0+LAMBDA1)/powf(Z1,.33333);
    ZETA0=Q+(1./(2.*VFERMI*VFERMI))*(1.-Q)*logf(1+powf((4.*L*VFERMI/1.919),2));
    if(E>1)tempf=logf(E);else tempf=0;
    ZETA=ZETA0*(1.+(1./(Z1*Z1))*(.08+.0015*Z2)*expf(-powf((7.6-tempf),2)));
    /*Calculate stopping for low YR*/
    tempf=VRMIN/powf(Z1,.6667); if(tempf<YRMIN)tempf=YRMIN;
    if(YR<=tempf){
      tempf=YRMIN*powf(Z1,.6667);if(VRMIN<tempf)VRMIN=tempf; /*VRMIN must be >= YRMIN^(2/3)*/
      tempf=VRMIN*VRMIN-0.8*VFERMI*VFERMI; if(tempf<0)tempf=0;
      VMIN=.5*(VRMIN+sqrtf(tempf));
      EEE=25*VMIN*VMIN;
      Sproton=pstop(Z2,EEE);
      /*Add Fermi Velocity Correction to Low Energy value*/
      if(EEE<=9999)EION=EEE;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
      for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
      if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
      VFCORR0=scoef[Z2][j];
      VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
      /*
      Following corrects for low-energy stopping, where little data exists.
      Traditionally, this is velocity-proportional stopping, however for
      light ions, light targets and semiconductors, a large variation exists.
      Note: HIPOWER down = Se up �� LAMBDA down = Se down
      */
      if(Z1==3)HIPOWER= 0.55;
      else if(Z2<7)HIPOWER= 0.375;
      else if(Z1<18 && (Z2==14 || Z2==32))HIPOWER= 0.375;
      else HIPOWER=.47; /*DEFAULT POWER IF NO SPECIAL CASE. TRIM-88 used 0.5*/
      SE[i]=Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1)*powf((E/EEE),HIPOWER);
      }
    /*Calculate stopping for NOT low YR*/
    else{
      Sproton=pstop(Z2,E);
      /*Add Fermi Velocity Correction - 1995*/
      if(E<=9999)EION=E;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
      for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
      if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
      VFCORR0=scoef[Z2][j];
      VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
      SE[i]=Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1);
      }
 /*
 */

    }
  /*SE(I)=SE(I)*10     'This converts Stopping in eV/(1E15-cm2) to eV-A2  */
  }
//for(i=0;i<nEner;i++)printf("---  %e\t%e\t%e\n",Ener[i],SE[i],pstop(14,Ener[i]));
}


void stop96d(int Z1, int Z2,double Ionweight, double *Ener, int nEner, double *SE){
/*Note this is only a version of  stop96, for programs using
double precission instead of "float". The internal calcs are still
 done in single precission!! */
/*
Z1=incident ion atomic number
Z2=target atomic number
M1=Weight of incident ion (if M1<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
nEner=dim of Ener as well as of SE
SE=pointer to a vector of nEner components where the resultin stoppings will be stored (ev/1e15atcm-2)
*/

float tempf,VFERMI,E,HE0,HE,B,HEH,YRMIN,VRMIN,V,VR,YR,Q,LAMBDA0,LAMBDA1,L,ZETA0,ZETA,VMIN,EEE,Sproton,EION,VFCORR0,VFCORR1,HIPOWER,M1;
int i,j;

/*Do checks*/
if(Z1<1 || Z1>92) {fprintf(stderr, "\nError: Z1 must be between 1 and 92 (%d)\n",Z1);exit(1);}
if(Z2<1 || Z2>92) {fprintf(stderr, "\nError: Z2 must be between 1 and 92 (%d)\n",Z2);exit(1);}
if(nEner<1){fprintf(stderr, "\nError: I need at least one energy  \n");exit(1);}

/*Fill some variables with data*/
VFERMI= scoef[Z2][7];       /*VFermi of Target*/
if(Ionweight<1.)M1=scoef[Z1][3] ;else M1=(float)Ionweight;    /*MOI weight if M1 wasn't specified*/


/*Loop for stopping values*/
for(i=0;i<nEner;i++){
  E=(float)(Ener[i]/M1);

  if(Z1==1){ /*stop_H*/
    SE[i]=(double)pstop(Z2,E);
    }


  else if(Z1==2){ /*stop_He */
    HE0=1;                            /*Minimum limit for "medium" energy range for He (units=keV/amu)*/
    if(HE0>E)HE=HE0;else HE=E;        /*HE=max{E,HE0}*/
    B=logf(HE);
    tempf=.2865+.1266*B-.001429*B*B+.02402*B*B*B-.01135*powf(B,4)+.001475*powf(B,5);
    if(tempf>30) tempf=30;
    HEH=1.-expf(-tempf);
    /* ADD Z1^3 EFFECT TO HE/H STOPPING POWER RATIO, HEH */
    if(HE<1)HE=1;  /*I think this assertion is innecessary!!!! (it its in the original code but seems redundant)???????*/
    tempf=(1.+(.007+.00005*Z2)*expf(-powf(7.6-logf(HE),2)));
    HEH*=(tempf*tempf);
    SE[i]=(double)(pstop(Z2,HE)*HEH*4.);
    /* If E was lower than HE0 , CALC. HE VELOCITY PROPORTIONAL STOPPING (low energy range)*/
    if(HE0>E)SE[i]*=(double)sqrtf(E/HE0);
    }


  else if(Z1>2){ /*stop_HI*/
    /*USE VELOCITY STOPPING FOR (YRMIN=VR/Z1**.67) <= 0.13  _OR_  VR<=1.0*/
    YRMIN=0.13;   /*Minimum limit for YR*/
    VRMIN=1.0;    /*Minimum limit for VR*/
    V=sqrtf(E/25.)/VFERMI;  /*V is the "relative velocity" but...relative to what????? (Fermi, Bohr,...????)*/
    if(V<1) VR=(3.*VFERMI/4.)*(1+(2.*V*V/3.)-powf(V,4.)/15.);
    else VR=V*VFERMI*(1.+1./(5.*V*V));
    /*put a minimum limit to VR*/
    if(VR<VRMIN)VR=VRMIN;
    /*calc YR as VR/Z1^(2/3) and put a minimum limit to it*/
    YR=VR/powf(Z1,.6667);
    if(YR<YRMIN)YR=YRMIN;
    /*Calculate the Ionization level (Q) of the ion at velocity YR*/
    tempf=-.803*powf(YR,0.3)+1.3167*powf(YR,0.6)+.38157*YR+.008983*YR*YR;
    if(tempf>50.)Q=1; else if(tempf<=0)Q=0; else  Q=1.-expf(-tempf);    /*It prevents underflow and wasting of time  in the Q calc. */
    /*NOW WE CONVERT IONIZATION LEVEL TO EFFECTIVE CHARGE
    ---------------------- Screening Distance of Ion (Lambda in B.& K.)*/
    for(j=22; j<=39 && scoef[93][j]<Q ;j++);j-=1;  /*Find Q Interpolation in SCOEF(93,22-39)*/
    if(j<22)j=22;else if(j>38)j=38; /*Boilerplate ??*/
    LAMBDA0=scoef[Z1][j];
    LAMBDA1=(Q-scoef[93][j])*(scoef[Z1][j+1]-scoef[Z1][j])/(scoef[93][j+1]-scoef[93][j]);
    L=(LAMBDA0+LAMBDA1)/powf(Z1,.33333);
    ZETA0=Q+(1./(2.*VFERMI*VFERMI))*(1.-Q)*logf(1+powf((4.*L*VFERMI/1.919),2));
    if(E>1)tempf=logf(E);else tempf=0;
    ZETA=ZETA0*(1.+(1./(Z1*Z1))*(.08+.0015*Z2)*expf(-powf((7.6-tempf),2)));
    /*Calculate stopping for low YR*/
    tempf=VRMIN/powf(Z1,.6667); if(tempf<YRMIN)tempf=YRMIN;
    if(YR<=tempf){
      tempf=YRMIN*powf(Z1,.6667);if(VRMIN<tempf)VRMIN=tempf; /*VRMIN must be >= YRMIN^(2/3)*/
      tempf=VRMIN*VRMIN-0.8*VFERMI*VFERMI; if(tempf<0)tempf=0;
      VMIN=.5*(VRMIN+sqrtf(tempf));
      EEE=25*VMIN*VMIN;
      Sproton=pstop(Z2,EEE);
      /*Add Fermi Velocity Correction to Low Energy value*/
      if(EEE<=9999)EION=EEE;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
      for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
      if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
      VFCORR0=scoef[Z2][j];
      VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
      /*
      Following corrects for low-energy stopping, where little data exists.
      Traditionally, this is velocity-proportional stopping, however for
      light ions, light targets and semiconductors, a large variation exists.
      Note: HIPOWER down = Se up �� LAMBDA down = Se down
      */
      if(Z1==3)HIPOWER= 0.55;
      else if(Z2<7)HIPOWER= 0.375;
      else if(Z1<18 && (Z2==14 || Z2==32))HIPOWER= 0.375;
      else HIPOWER=.47; /*DEFAULT POWER IF NO SPECIAL CASE. TRIM-88 used 0.5*/
      SE[i]=(double)(Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1)*powf((E/EEE),HIPOWER));
      }
    /*Calculate stopping for NOT low YR*/
    else{
      Sproton=pstop(Z2,E);
      /*Add Fermi Velocity Correction - 1995*/
      if(E<=9999)EION=E;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
      for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
      if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
      VFCORR0=scoef[Z2][j];
      VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
      SE[i]=(double)(Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1));
      }
 /*
 */

    }
  /*SE(I)=SE(I)*10     'This converts Stopping in eV/(1E15-cm2) to eV-A2  */
  }
//for(i=0;i<nEner;i++)printf("---  %e\t%e\t%e\n",Ener[i],SE[i],pstop(14,Ener[i]));
}

double nuclearstopping_ZBL(int Z1, int Z2,double Ionweight, double trgtweight, double *Ener, int nEner, double *SN)
{
/*
nuclearstopping_ZBL calculates ZBL 1990 nuclear stopping powers
*/
/*
Z1=incident ion atomic number
Z2=target atomic number
Ionweight=Weight of incident ion (if Ionweight<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
trgtweight=Weight of target atom (if trgtweight<1, the Most Abundant Isotope weight associated to Z2 in SCOEF.95A will be used)
Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
nEner=dim of Ener as well as of SE
SN=pointer to a vector of nEner components where the resulting stoppings will be stored (ev/1e15atcm-2)

Return value: spmax: maximum nuclear stopping power computed. If nEner=1, then spmax=SN[0] !!!

IMPORTANT: SN and SP[] are returned in eV/(1e15at/cm2). Be careful of converting to keV/(1e15at/cm2) when using the value for internal Hotstop calcs
*/
double M1,M2,EPSIL,temp,spmax;
int i;

/*Determine ion and target massess*/
if(trgtweight<1.)M2=(double)scoef[Z2][3] ;else M2=(double)trgtweight;    /*MOI weight if M2 wasn't specified*/
if(Ionweight<1.)M1=(double)scoef[Z1][3] ;else M1=(double)Ionweight;    /*MOI weight if M1 wasn't specified*/

for(i=0,spmax=0;i<nEner;i++){
  /*Epsilon is the reduced energy of the ion/target combination.*/
  temp=(M1+M2)*(pow(Z1,0.23)+pow(Z2,0.23));
  EPSIL=32.53*M2*Ener[i]/(Z1*Z2*temp);
  if(EPSIL==0)SN[i]=0;
  else if(EPSIL>=30) SN[i]=log(EPSIL)/(2*EPSIL);
  else SN[i]=.5*log(1+1.1383*EPSIL)/(EPSIL+(.01321*pow(EPSIL,.21226))+(.19593*sqrt(EPSIL)));
  /*convert from LSS reduced units to eV-cm2/1E15 */
  SN[i]*=Z1*Z2*M1*8.462/temp;
  /*Find the maximum SN*/
  if(SN[i]>spmax)spmax=SN[i];
  }
return(spmax);
}

/*
int main(void)
{
float E,Sproton;
int i,j,NE;
float *S,*Ener;
readstopcoef();
NE=1000;
if(!(S=calloc(NE,sizeof(float)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
if(!(Ener=calloc(NE,sizeof(float)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
for(i=0,E=0;i<NE;i++,E+=1.)Ener[i]=E;
stop96(1,16,Ener,NE,S);
for(i=0;i<NE;i++)printf("%e\t%e\n",Ener[i],S[i]);
}
*/

---