Main Page | Alphabetical List | Data Structures | Directories | File List | Data Fields | Globals | Related Pages

stop96.c

Go to the documentation of this file.
00001 /*
00002     STOP96.c: Ion stopping subroutines used by hotstop
00003     By Carlos Pascual-Izarra
00004 
00005  Copyright (C) 2002 Carlos Pascual-Izarra
00006 
00007    This program is free software; you can redistribute it and/or modify
00008    it under the terms of the GNU General Public License as published by
00009    the Free Software Foundation; either version 2, or (at your option)
00010    any later version.
00011 
00012    This program is distributed in the hope that it will be useful,
00013    but WITHOUT ANY WARRANTY; without even the implied warranty of
00014    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015    GNU General Public License for more details.
00016 
00017    You should have received a copy of the GNU General Public License
00018    along with this program; if not, write to the Free Software Foundation,
00019    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00020 
00021     This is a C version of the STOP96 subroutine by J.F.Ziegler for 
00022     being used in TRIM / SRIM. The original from Ziegler was written
00023     in BASIC.
00024 
00025     IMPORTANT:
00026     It must be clear that J.F.Ziegler or any other person related to 
00027     TRIM/SRIM have absolutely no responsability of any limitation or
00028     bug related to the present code.
00029     
00030     The following is the license notice of SRIM:
00031     
00032     "SRIM.EXE, (C) 1984-2002, James F. Ziegler
00033     Permission to use, copy, modify and distribute 
00034     this software for any non-commercial purpose and 
00035     without fee is hereby granted, provided that 
00036     this copyright and permission notice appear
00037     on all copies of the software. The name of 
00038     the author may not be used in any advertising 
00039     or publicity pertaining to the use of the software.
00040     The author makes no warranty or representations 
00041     about the suitability of the software for any 
00042     purpose. It is provided "AS IS" without any express 
00043     or implied warranty, including the implied warranties 
00044     of merchantability, fitness for a particular purpose 
00045     and non-infringement. The author shall not be liable 
00046     for any direct, indirect, special or consequential 
00047     damages resulting from the loss of use, data or 
00048     projects, whether in an action of contract or tort, 
00049     arising out of or in connection with the use or 
00050     performance of this software. Downloading any 
00051     software from this site is implicit acceptance of 
00052     these conditions. The term "SRIM" as used to 
00053     refer to topics concerning the  Stopping and 
00054     Range of Ions in Matter is copyright registered 
00055     in the United States Copyright Office and may 
00056     not be used for commercial purposes without permission. 
00057     For commercial applications, licenses may be obtained 
00058     for incorporating SRIM into applications by writing to 
00059     SRIM.org."  
00060 
00061   
00062 
00063  */
00064  /**\file stop96.c
00065  stop96.c contains functions related to stopping force calculations. This file has been taken from the Hotstop source code ( http://hotstop.sourceforge.net ).
00066  
00067 */
00068 #include <stdio.h>
00069 #include <stdlib.h>
00070 #include <math.h>     /* REMEMBER: use -lm when compiling */
00071 #include "stop96.h"  /*Contains the definitions for the stop96.c functions*/
00072 #include "compilopt.h"
00073 
00074 
00075 void readstopcoef(char *afilename, char *bfilename)
00076 {
00077 /*Fills the scoef[94,55] array with the data in SCOEF.95A and SCOEF.95B files*/
00078 /* scoef[94,55] must be previously declared as a global variable */
00079 /*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)*/
00080 
00081 FILE *afile, *bfile;
00082 int i,j;
00083 char stemp[300];
00084 
00085 /*Open files*/
00086 if (!(afile=fopen(afilename, "r"))){fprintf(stderr, "\ncannot open '%s' for read\n",afilename);exit(1);}
00087 if (!(bfile=fopen(bfilename, "r"))){fprintf(stderr, "\ncannot open '%s' for read\n",bfilename);exit(1);}
00088 
00089 /*Skip headers*/
00090 for(i=0;i<2;i++) fgets(stemp,300,afile) ;  /*Skip Headers in file A*/
00091 for(i=0;i<2;i++) fgets(stemp,300,bfile) ;  /*Skip Headers in file B*/
00092 
00093 /*fill scoef*/
00094 for(i=1;i<=93;i++) for(j=1;j<=16;j++)fscanf(afile,"%f",&scoef[i][j]); 
00095 for(i=1;i<=93;i++) for(j=17;j<=54;j++)fscanf(bfile,"%f",&scoef[i][j]);
00096 
00097 
00098 /*Close files*/
00099 fclose(afile);
00100 fclose(bfile);
00101 }
00102 
00103 float pstop(int Z2, float E)
00104 {
00105 
00106 /* scoef[94,55] must be previously declared as a global variable */
00107 float X,SL,SH,ppower,Sproton;
00108 float PE0;
00109 PE0=10;
00110 /*High (E>10000keV/amu) energy stopping*/
00111 if(E>1e4){
00112   X=logf(E)/E;
00113   Sproton=scoef[Z2][17]+(scoef[Z2][18]*X)+(scoef[Z2][19]*X*X)+(scoef[Z2][20]/X);
00114   }
00115 /*Medium energies: (PE0<E<1e4)  , PE0=is the low limit for "medium" energy range (in kev/amu)*/
00116 else if (E>PE0){
00117   SL=scoef[Z2][9]*powf(E,scoef[Z2][10])+scoef[Z2][11]*powf(E,scoef[Z2][12]);
00118   SH=scoef[Z2][13]/powf(E,scoef[Z2][14])*logf((scoef[Z2][15]/E)+scoef[Z2][16]*E);
00119   Sproton=SL*SH/(SL+SH);
00120   }
00121 /*Low energies (E<PE0) : approx velocity proportional stopping power, conveniently scaled for continuity with medium E range*/  
00122 else if(E>0){
00123   if(Z2<7)ppower=0.35; else ppower=0.45;
00124   SL=scoef[Z2][9]*powf(PE0,scoef[Z2][10])+scoef[Z2][11]*powf(PE0,scoef[Z2][12]);
00125   SH=scoef[Z2][13]/powf(PE0,scoef[Z2][14])*logf((scoef[Z2][15]/PE0)+scoef[Z2][16]*PE0);
00126   Sproton=(SL*SH/(SL+SH))*powf(E/PE0,ppower);
00127   }
00128 else Sproton=0.;
00129 
00130 return(Sproton);
00131 }
00132 
00133 
00134 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)
00135 {
00136 /*This function is identical to pstop but it gets all the parameters from its arguments instead of looking at scoef[Z2][x]*/
00137 float X,SL,SH,ppower,Sproton;
00138 float PE0;
00139 PE0=10;
00140 /*High (E>10000keV/amu) energy stopping*/
00141 if(E>1e4){
00142   X=logf(E)/E;
00143   Sproton=HE0+(HE1*X)+(HE2*X*X)+(HE3/X);
00144   }
00145 /*Medium energies: (PE0<E<1e4)  , PE0=is the low limit for "medium" energy range (in kev/amu)*/
00146 else if (E>PE0){
00147   SL=P0*powf(E,P1)+P2*powf(E,P3);
00148   SH=P4/powf(E,P5)*logf((P6/E)+P7*E);
00149   Sproton=SL*SH/(SL+SH);
00150   }
00151 /*Low energies (E<PE0) : approx velocity proportional stopping power, conveniently scaled for continuity with medium E range*/  
00152 else if(E>0){
00153   if(Z2<7)ppower=0.35; else ppower=0.45;
00154   SL=P0*powf(PE0,P1)+P2*powf(PE0,P3);
00155   SH=P4/powf(PE0,P5)*logf((P6/PE0)+P7*PE0);
00156   Sproton=(SL*SH/(SL+SH))*powf(E/PE0,ppower);
00157   }
00158 else Sproton=0.;
00159 
00160 return(Sproton);
00161 }
00162 
00163 void stop96(int Z1, int Z2,float Ionweight, float *Ener, int nEner, float *SE){
00164 
00165 /*
00166 Z1=incident ion atomic number
00167 Z2=target atomic number
00168 M1=Weight of incident ion (if M1<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
00169 Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
00170 nEner=dim of Ener as well as of SE
00171 SE=pointer to a vector of nEner components where the resultin stoppings will be stored (ev/1e15atcm-2)
00172 
00173 
00174 */
00175 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;
00176 int i,j;
00177 
00178 /*Do checks*/
00179 if(Z1<1 || Z1>92) {fprintf(stderr, "\nError: Z1 must be between 1 and 92 \n");exit(1);}
00180 if(Z2<1 || Z2>92) {fprintf(stderr, "\nError: Z2 must be between 1 and 92 \n");exit(1);}
00181 if(nEner<1){fprintf(stderr, "\nError: I need at least one energy  \n");exit(1);}
00182 
00183 
00184 /*Fill some variables with data*/
00185 VFERMI= scoef[Z2][7];       /*VFermi of Target*/
00186 if(Ionweight<1)M1=scoef[Z1][3]; else M1=Ionweight;    /*MOI weight if M1 wasn't specified*/
00187 
00188 /*Loop for stopping values*/
00189 for(i=0;i<nEner;i++){
00190   E=Ener[i]/M1;
00191 
00192   if(Z1==1){ /*stop_H*/
00193     SE[i]=pstop(Z2,E);
00194     }
00195 
00196 
00197   else if(Z1==2){ /*stop_He */
00198     HE0=1;                            /*Minimum limit for "medium" energy range for He (units=keV/amu)*/
00199     if(HE0>E)HE=HE0;else HE=E;        /*HE=max{E,HE0}*/
00200     B=logf(HE);
00201     tempf=.2865+.1266*B-.001429*B*B+.02402*B*B*B-.01135*powf(B,4)+.001475*powf(B,5);
00202     if(tempf>30) tempf=30;
00203     HEH=1.-expf(-tempf);
00204     /* ADD Z1^3 EFFECT TO HE/H STOPPING POWER RATIO, HEH */
00205     if(HE<1)HE=1;  /*I think this assertion is innecessary!!!! (it its in the original code but seems redundant)???????*/
00206     tempf=(1.+(.007+.00005*Z2)*expf(-powf(7.6-logf(HE),2)));
00207     HEH*=(tempf*tempf);
00208     SE[i]=pstop(Z2,HE)*HEH*4.;
00209     /* If E was lower than HE0 , CALC. HE VELOCITY PROPORTIONAL STOPPING (low energy range)*/
00210     if(HE0>E)SE[i]*=sqrtf(E/HE0);
00211     }
00212       
00213 
00214   else if(Z1>2){ /*stop_HI*/
00215     /*USE VELOCITY STOPPING FOR (YRMIN=VR/Z1**.67) <= 0.13  _OR_  VR<=1.0*/
00216     YRMIN=0.13;   /*Minimum limit for YR*/
00217     VRMIN=1.0;    /*Minimum limit for VR*/
00218     V=sqrtf(E/25.)/VFERMI;  /*V is the "relative velocity" but...relative to what????? (Fermi, Bohr,...????)*/
00219     if(V<1) VR=(3.*VFERMI/4.)*(1+(2.*V*V/3.)-powf(V,4.)/15.);
00220     else VR=V*VFERMI*(1.+1./(5.*V*V));
00221     /*put a minimum limit to VR*/
00222     if(VR<VRMIN)VR=VRMIN;
00223     /*calc YR as VR/Z1^(2/3) and put a minimum limit to it*/
00224     YR=VR/powf(Z1,.6667);
00225     if(YR<YRMIN)YR=YRMIN;
00226     /*Calculate the Ionization level (Q) of the ion at velocity YR*/
00227     tempf=-.803*powf(YR,0.3)+1.3167*powf(YR,0.6)+.38157*YR+.008983*YR*YR;
00228     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. */
00229     /*NOW WE CONVERT IONIZATION LEVEL TO EFFECTIVE CHARGE
00230     ---------------------- Screening Distance of Ion (Lambda in B.& K.)*/
00231     for(j=22; j<=39 && scoef[93][j]<Q ;j++);j-=1;  /*Find Q Interpolation in SCOEF(93,22-39)*/
00232     if(j<22)j=22;else if(j>38)j=38; /*Boilerplate ??*/
00233     LAMBDA0=scoef[Z1][j];
00234     LAMBDA1=(Q-scoef[93][j])*(scoef[Z1][j+1]-scoef[Z1][j])/(scoef[93][j+1]-scoef[93][j]);
00235     L=(LAMBDA0+LAMBDA1)/powf(Z1,.33333);
00236     ZETA0=Q+(1./(2.*VFERMI*VFERMI))*(1.-Q)*logf(1+powf((4.*L*VFERMI/1.919),2));
00237     if(E>1)tempf=logf(E);else tempf=0;
00238     ZETA=ZETA0*(1.+(1./(Z1*Z1))*(.08+.0015*Z2)*expf(-powf((7.6-tempf),2)));
00239     /*Calculate stopping for low YR*/
00240     tempf=VRMIN/powf(Z1,.6667); if(tempf<YRMIN)tempf=YRMIN;
00241     if(YR<=tempf){
00242       tempf=YRMIN*powf(Z1,.6667);if(VRMIN<tempf)VRMIN=tempf; /*VRMIN must be >= YRMIN^(2/3)*/
00243       tempf=VRMIN*VRMIN-0.8*VFERMI*VFERMI; if(tempf<0)tempf=0;
00244       VMIN=.5*(VRMIN+sqrtf(tempf));
00245       EEE=25*VMIN*VMIN;
00246       Sproton=pstop(Z2,EEE);
00247       /*Add Fermi Velocity Correction to Low Energy value*/
00248       if(EEE<=9999)EION=EEE;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
00249       for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
00250       if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
00251       VFCORR0=scoef[Z2][j];
00252       VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
00253       /*
00254       Following corrects for low-energy stopping, where little data exists.
00255       Traditionally, this is velocity-proportional stopping, however for
00256       light ions, light targets and semiconductors, a large variation exists.
00257       Note: HIPOWER down = Se up �� LAMBDA down = Se down 
00258       */
00259       if(Z1==3)HIPOWER= 0.55;
00260       else if(Z2<7)HIPOWER= 0.375;
00261       else if(Z1<18 && (Z2==14 || Z2==32))HIPOWER= 0.375;
00262       else HIPOWER=.47; /*DEFAULT POWER IF NO SPECIAL CASE. TRIM-88 used 0.5*/
00263       SE[i]=Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1)*powf((E/EEE),HIPOWER);
00264       }
00265     /*Calculate stopping for NOT low YR*/ 
00266     else{
00267       Sproton=pstop(Z2,E);
00268       /*Add Fermi Velocity Correction - 1995*/
00269       if(E<=9999)EION=E;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
00270       for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
00271       if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
00272       VFCORR0=scoef[Z2][j];
00273       VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
00274       SE[i]=Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1); 
00275       }
00276  /*  
00277  */
00278   
00279     }
00280   /*SE(I)=SE(I)*10     'This converts Stopping in eV/(1E15-cm2) to eV-A2  */
00281   }
00282 //for(i=0;i<nEner;i++)printf("---  %e\t%e\t%e\n",Ener[i],SE[i],pstop(14,Ener[i]));
00283 }
00284 
00285 
00286 void stop96d(int Z1, int Z2,double Ionweight, double *Ener, int nEner, double *SE){
00287 /*Note this is only a version of  stop96, for programs using 
00288 double precission instead of "float". The internal calcs are still
00289  done in single precission!! */
00290 /*
00291 Z1=incident ion atomic number
00292 Z2=target atomic number
00293 M1=Weight of incident ion (if M1<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
00294 Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
00295 nEner=dim of Ener as well as of SE
00296 SE=pointer to a vector of nEner components where the resultin stoppings will be stored (ev/1e15atcm-2)
00297 */
00298 
00299 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;
00300 int i,j;
00301 
00302 /*Do checks*/
00303 if(Z1<1 || Z1>92) {fprintf(stderr, "\nError: Z1 must be between 1 and 92 (%d)\n",Z1);exit(1);}
00304 if(Z2<1 || Z2>92) {fprintf(stderr, "\nError: Z2 must be between 1 and 92 (%d)\n",Z2);exit(1);}
00305 if(nEner<1){fprintf(stderr, "\nError: I need at least one energy  \n");exit(1);}
00306 
00307 /*Fill some variables with data*/
00308 VFERMI= scoef[Z2][7];       /*VFermi of Target*/
00309 if(Ionweight<1.)M1=scoef[Z1][3] ;else M1=(float)Ionweight;    /*MOI weight if M1 wasn't specified*/
00310 
00311 
00312 /*Loop for stopping values*/
00313 for(i=0;i<nEner;i++){
00314   E=(float)(Ener[i]/M1);
00315 
00316   if(Z1==1){ /*stop_H*/
00317     SE[i]=(double)pstop(Z2,E);
00318     }
00319 
00320 
00321   else if(Z1==2){ /*stop_He */
00322     HE0=1;                            /*Minimum limit for "medium" energy range for He (units=keV/amu)*/
00323     if(HE0>E)HE=HE0;else HE=E;        /*HE=max{E,HE0}*/
00324     B=logf(HE);
00325     tempf=.2865+.1266*B-.001429*B*B+.02402*B*B*B-.01135*powf(B,4)+.001475*powf(B,5);
00326     if(tempf>30) tempf=30;
00327     HEH=1.-expf(-tempf);
00328     /* ADD Z1^3 EFFECT TO HE/H STOPPING POWER RATIO, HEH */
00329     if(HE<1)HE=1;  /*I think this assertion is innecessary!!!! (it its in the original code but seems redundant)???????*/
00330     tempf=(1.+(.007+.00005*Z2)*expf(-powf(7.6-logf(HE),2)));
00331     HEH*=(tempf*tempf);
00332     SE[i]=(double)(pstop(Z2,HE)*HEH*4.);
00333     /* If E was lower than HE0 , CALC. HE VELOCITY PROPORTIONAL STOPPING (low energy range)*/
00334     if(HE0>E)SE[i]*=(double)sqrtf(E/HE0);
00335     }
00336       
00337 
00338   else if(Z1>2){ /*stop_HI*/
00339     /*USE VELOCITY STOPPING FOR (YRMIN=VR/Z1**.67) <= 0.13  _OR_  VR<=1.0*/
00340     YRMIN=0.13;   /*Minimum limit for YR*/
00341     VRMIN=1.0;    /*Minimum limit for VR*/
00342     V=sqrtf(E/25.)/VFERMI;  /*V is the "relative velocity" but...relative to what????? (Fermi, Bohr,...????)*/
00343     if(V<1) VR=(3.*VFERMI/4.)*(1+(2.*V*V/3.)-powf(V,4.)/15.);
00344     else VR=V*VFERMI*(1.+1./(5.*V*V));
00345     /*put a minimum limit to VR*/
00346     if(VR<VRMIN)VR=VRMIN;
00347     /*calc YR as VR/Z1^(2/3) and put a minimum limit to it*/
00348     YR=VR/powf(Z1,.6667);
00349     if(YR<YRMIN)YR=YRMIN;
00350     /*Calculate the Ionization level (Q) of the ion at velocity YR*/
00351     tempf=-.803*powf(YR,0.3)+1.3167*powf(YR,0.6)+.38157*YR+.008983*YR*YR;
00352     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. */
00353     /*NOW WE CONVERT IONIZATION LEVEL TO EFFECTIVE CHARGE
00354     ---------------------- Screening Distance of Ion (Lambda in B.& K.)*/
00355     for(j=22; j<=39 && scoef[93][j]<Q ;j++);j-=1;  /*Find Q Interpolation in SCOEF(93,22-39)*/
00356     if(j<22)j=22;else if(j>38)j=38; /*Boilerplate ??*/
00357     LAMBDA0=scoef[Z1][j];
00358     LAMBDA1=(Q-scoef[93][j])*(scoef[Z1][j+1]-scoef[Z1][j])/(scoef[93][j+1]-scoef[93][j]);
00359     L=(LAMBDA0+LAMBDA1)/powf(Z1,.33333);
00360     ZETA0=Q+(1./(2.*VFERMI*VFERMI))*(1.-Q)*logf(1+powf((4.*L*VFERMI/1.919),2));
00361     if(E>1)tempf=logf(E);else tempf=0;
00362     ZETA=ZETA0*(1.+(1./(Z1*Z1))*(.08+.0015*Z2)*expf(-powf((7.6-tempf),2)));
00363     /*Calculate stopping for low YR*/
00364     tempf=VRMIN/powf(Z1,.6667); if(tempf<YRMIN)tempf=YRMIN;
00365     if(YR<=tempf){
00366       tempf=YRMIN*powf(Z1,.6667);if(VRMIN<tempf)VRMIN=tempf; /*VRMIN must be >= YRMIN^(2/3)*/
00367       tempf=VRMIN*VRMIN-0.8*VFERMI*VFERMI; if(tempf<0)tempf=0;
00368       VMIN=.5*(VRMIN+sqrtf(tempf));
00369       EEE=25*VMIN*VMIN;
00370       Sproton=pstop(Z2,EEE);
00371       /*Add Fermi Velocity Correction to Low Energy value*/
00372       if(EEE<=9999)EION=EEE;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
00373       for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
00374       if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
00375       VFCORR0=scoef[Z2][j];
00376       VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
00377       /*
00378       Following corrects for low-energy stopping, where little data exists.
00379       Traditionally, this is velocity-proportional stopping, however for
00380       light ions, light targets and semiconductors, a large variation exists.
00381       Note: HIPOWER down = Se up �� LAMBDA down = Se down 
00382       */
00383       if(Z1==3)HIPOWER= 0.55;
00384       else if(Z2<7)HIPOWER= 0.375;
00385       else if(Z1<18 && (Z2==14 || Z2==32))HIPOWER= 0.375;
00386       else HIPOWER=.47; /*DEFAULT POWER IF NO SPECIAL CASE. TRIM-88 used 0.5*/
00387       SE[i]=(double)(Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1)*powf((E/EEE),HIPOWER));
00388       }
00389     /*Calculate stopping for NOT low YR*/ 
00390     else{
00391       Sproton=pstop(Z2,E);
00392       /*Add Fermi Velocity Correction - 1995*/
00393       if(E<=9999)EION=E;else EION=9999; /*Correction is only valid for E <1E4 keV/amu*/
00394       for(j=41;j<=53 && scoef[93][j]<EION ;j++);j-=1;  /*Find E Interpolation in SCOEF(93,41-54)*/
00395       if(j<41)j=41;else if(j>53)j=53; /*Boilerplate ??*/
00396       VFCORR0=scoef[Z2][j];
00397       VFCORR1=(EION-scoef[93][j])*(scoef[Z2][j+1]-scoef[Z2][j])/(scoef[93][j+1]-scoef[93][j]);
00398       SE[i]=(double)(Sproton*powf(ZETA*Z1,2)*(VFCORR0+VFCORR1)); 
00399       }
00400  /*  
00401  */
00402   
00403     }
00404   /*SE(I)=SE(I)*10     'This converts Stopping in eV/(1E15-cm2) to eV-A2  */
00405   }
00406 //for(i=0;i<nEner;i++)printf("---  %e\t%e\t%e\n",Ener[i],SE[i],pstop(14,Ener[i]));
00407 }
00408 
00409 double nuclearstopping_ZBL(int Z1, int Z2,double Ionweight, double trgtweight, double *Ener, int nEner, double *SN)
00410 {
00411 /*
00412 nuclearstopping_ZBL calculates ZBL 1990 nuclear stopping powers
00413 */
00414 /*
00415 Z1=incident ion atomic number
00416 Z2=target atomic number
00417 Ionweight=Weight of incident ion (if Ionweight<1, the Most Abundant Isotope weight associated to Z1 in SCOEF.95A will be used)
00418 trgtweight=Weight of target atom (if trgtweight<1, the Most Abundant Isotope weight associated to Z2 in SCOEF.95A will be used)
00419 Ener= pointer to vector of nEner components storing the list of energies wanted (in keV)
00420 nEner=dim of Ener as well as of SE
00421 SN=pointer to a vector of nEner components where the resulting stoppings will be stored (ev/1e15atcm-2)
00422 
00423 Return value: spmax: maximum nuclear stopping power computed. If nEner=1, then spmax=SN[0] !!!
00424 
00425 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
00426 */
00427 double M1,M2,EPSIL,temp,spmax;
00428 int i;
00429 
00430 /*Determine ion and target massess*/
00431 if(trgtweight<1.)M2=(double)scoef[Z2][3] ;else M2=(double)trgtweight;    /*MOI weight if M2 wasn't specified*/
00432 if(Ionweight<1.)M1=(double)scoef[Z1][3] ;else M1=(double)Ionweight;    /*MOI weight if M1 wasn't specified*/
00433 
00434 for(i=0,spmax=0;i<nEner;i++){
00435   /*Epsilon is the reduced energy of the ion/target combination.*/
00436   temp=(M1+M2)*(pow(Z1,0.23)+pow(Z2,0.23));
00437   EPSIL=32.53*M2*Ener[i]/(Z1*Z2*temp);
00438   if(EPSIL==0)SN[i]=0;
00439   else if(EPSIL>=30) SN[i]=log(EPSIL)/(2*EPSIL); 
00440   else SN[i]=.5*log(1+1.1383*EPSIL)/(EPSIL+(.01321*pow(EPSIL,.21226))+(.19593*sqrt(EPSIL)));
00441   /*convert from LSS reduced units to eV-cm2/1E15 */
00442   SN[i]*=Z1*Z2*M1*8.462/temp;
00443   /*Find the maximum SN*/
00444   if(SN[i]>spmax)spmax=SN[i];
00445   }
00446 return(spmax);
00447 }
00448 
00449 /*
00450 int main(void)
00451 {
00452 float E,Sproton;
00453 int i,j,NE;
00454 float *S,*Ener;
00455 readstopcoef(); 
00456 NE=1000;
00457 if(!(S=calloc(NE,sizeof(float)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
00458 if(!(Ener=calloc(NE,sizeof(float)))){fprintf(stderr,"\n Error allocating memory\n");exit(1);}
00459 for(i=0,E=0;i<NE;i++,E+=1.)Ener[i]=E;
00460 stop96(1,16,Ener,NE,S);
00461 for(i=0;i<NE;i++)printf("%e\t%e\n",Ener[i],S[i]);
00462 }
00463 */

Generated on Fri Jul 15 20:43:49 2005 for LibCPIXE API by  doxygen 1.3.9.1