LMU/TRIMSP
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
d76d40ea4ce638ff8cf05012a2f2bf0fb7ad48e4
TRIMSP/fortran/trimspNL.F
T

5572 lines
202 KiB
FortranFixed
Raw Blame History

C Created June 2000
C
C Expanded to heterostructures with N Layer and MAXEL elements
C Zaher Salman <zaher.salman@psi.ch>
C
C Copyright W. Eckstein, IPP GARCHING, FRG
C Copyright Z. Salman and the LEM group PSI
C if you use this code please cite
C
C Eckstein, W. Computer Simulation of Ion-Solid Interactions;
C Springer Series in Materials Science;
C Springer-Verlag: Berlin Heidelberg, 1991.
C
C Historical overview:
C PROGRAM TRVMC95 VERSION AT IPP GARCHING NOVEMBER 1995
C MOMENTS OF DISTRIBUTIONS (RANGE, ENERGY AND
C ANGLE OF BACKSCATTERED AND SPUTTERED ATOMS)
C INTRODUCTION OF TAUPSI, TAUPSIR
C
C (PROGRAM TRVMC VERSION AT IPP GARCHING MARCH 1993)
C
C VECTORIZED TRIM FOR SPUTTERING, MULTI-COMPONENT TARGET
C 3 LAYERS A 5 COMPONENTS
C BACKWARD AND TRANSMISSION SPUTTERING
C
C W.ECKSTEIN IPP/PWW GARCHING CRAY-XMP
C WORKSTATIONS
C
C BASED ON TRSP1CN, TRIMSP3D (W.ECKSTEIN AND J.P.BIERSACK)
C VECTORIZATION BASED ON TRIM.VEC (M.BASKES SANDIA LIVERMORE)
C
C BACKSCATTERING AND TRANSMISSION OF PROJECTILES
C BACKWARD AND TRANSMISSION SPUTTERING
C ENERGY DISTRIBUTIONS IN SPECIFIC DIRECTIONS
C ENERGY DISTRIBUTIONS IN 100 STEPS UP TO INCIDENT ENERGY
C ANGULAR DISTRIBUTIONS IN STEPS OF COS = 0.05 (CONST. SOLID ANGLE)
C
c-------------------------------------------
c check OS
c
#if defined( _WIN32 )
#define OS_WIN
#else
#define OS_UNIX
#endif
IMPLICIT NONE
C=======================================================================
C DECLARATIONS AND GLOBAL STORAGE
C=======================================================================
CHARACTER*16 TRIMSP_VERSION
C These parameters are related to the maximum number of layers MAXNL
C and maximum number of points in the depth distribution MAXD
INTEGER MAXD,MAXD1,MAXD2,MAXD5,MAXDNL5,MAXD2p2,MAXD2meagb
& ,MAXD2meab,MAXEL
INTEGER MAXNL,MAXNL5,MAXNLp25,MAXNL5p2,MAXNLm15
C This is the only point where the number of layers and depth
C profile are changed. All other parameters should be changed
C accordingly.
PARAMETER (MAXD=500)
PARAMETER (MAXNL=60)
C Maximum number of elements in each layer, was limited to 5.
PARAMETER (MAXEL=12)
PARAMETER (MAXD1=MAXD+1)
PARAMETER (MAXD2=MAXD+2)
PARAMETER (MAXD5=MAXD*MAXEL)
PARAMETER (MAXD2p2=MAXD2*MAXD2)
PARAMETER (MAXD2meagb=MAXD2*36*22)
PARAMETER (MAXD2meab=MAXD2*22)
PARAMETER (MAXNL5=MAXNL*MAXEL)
PARAMETER (MAXNLp25=MAXNL*MAXNL5)
PARAMETER (MAXDNL5=MAXNL*MAXD5)
PARAMETER (MAXNL5p2=MAXNL5*MAXNL5*MAXD)
PARAMETER (MAXNLm15=(MAXNL-1)*MAXEL)
PARAMETER (TRIMSP_VERSION='1.3.2')
C
C Named integer options used by the old input format. The input
C file still contains the numeric values; these PARAMETER names make
C comparisons and output descriptions easier to read.
C
INTEGER POT_KRC,POT_MOLIERE,POT_ZBL
INTEGER STOP_LS,STOP_OR,STOP_MIXED,STOP_ICRU49,STOP_ZIEGLER
INTEGER IRL_OFF
PARAMETER (POT_KRC=1,POT_MOLIERE=2,POT_ZBL=3)
PARAMETER (STOP_LS=1,STOP_OR=2,STOP_MIXED=3)
PARAMETER (STOP_ICRU49=4,STOP_ZIEGLER=5)
PARAMETER (IRL_OFF=0)
LOGICAL TEST(64),TESTR(2000),TEST1(2000)
LOGICAL EQUAL
INTEGER*4 ISRCHFGT,ISRCHFGE,ILLZ
INTEGER N,L,LL,NH,NUM,KK
INTEGER I,J,IV
INTEGER tryE,negE
INTEGER depth_interval_flag
INTEGER*4 seconds_start_total,seconds_stop_total
INTEGER*4 NREC1,NREC2,NE1,K,NGIK,ICW
INTEGER*4 ISEED,ISEED2,ISEED3
INTEGER*4 JJR(2000,2),INOUT(2000,2),LRR(2000,2)
INTEGER*4 IDMAX(2000),IKR(2000)
INTEGER*4 number_in_layer(MAXNL),laufzahl
INTEGER*4 IRP(0:MAXD1),IPL(MAXD),IPLB(MAXD),IPLT(MAXD)
INTEGER*4 ICD(MAXD,MAXNL5),ICDT(MAXD),ICDJT(MAXNL5)
& ,ICDIRJ(MAXNL5,MAXNL5),ICDR(MAXD,MAXNL5),ICDTR(MAXD)
& ,ICDJTR(MAXNL5),ICDIRI(MAXD,MAXNL5,MAXNL5)
& ,ICDIRN(MAXD,MAXNL5),ICDITR(MAXNL5)
INTEGER*4 KADB(20),KADT(20),KADS(20),KADST(20) ,KADRIP(20,MAXNLm15
& ),KADRIS(20,MAXNLm15),KADROP(20,MAXNLm15),KADROS(20,MAXNLm15)
& ,KADSJ(20,MAXNLm15),KADSL(20,6),KDSTJ(20,MAXNLm15),KDSTL(20,6
& )
INTEGER*4 IBSP(MAXNL5),ISPAL(MAXNL),IBSPL(MAXNL5),ISPIP(MAXNL5)
& ,ISPIS(MAXNL5),ISPOP(MAXNL5),ISPOS(MAXNL5)
INTEGER*4 ITSP(MAXNL5),ISPALT(MAXNL)
& ,ISPIPT(MAXNL5),ISPIST(MAXNL5),ISPOPT(MAXNL5),ISPOST(MAXNL5)
INTEGER*4 KO(600,MAXNL5,2)
INTEGER*4 MEAB(MAXD2,22),MAGB(62,22),MEAGB(MAXD2,36,22)
& ,MEABL(75,21),MEPB(MAXD2,MAXD2)
INTEGER*4 MEAT(MAXD2,22),MAGT(62,22),MEAGT(MAXD2,36,22),
& MEATL(75,21),MEPT(MAXD2,MAXD2)
INTEGER*4 MEAS(MAXD2,22,MAXNLm15),MAGS(62,22,MAXNLm15),MAGSA(62,32
& ,MAXNLm15),MEAGS(MAXD2,12,22,MAXNLm15) ,MEASL(75,21,MAXNLm15)
INTEGER*4 MEAST(MAXD2,22,MAXNLm15),MAGST(62,22,MAXNLm15),MEASTL(75
& ,21,MAXNLm15)
C NJ(I) is the number of elements in layer I
INTEGER*4 NJ(MAXNL),JT(MAXNL),ILD(MAXNL),NJJ
INTEGER*4 LLL(64),JJJ(64),IK(64)
INTEGER*4 me(5000),nli(0:MAXNL),irpl(MAXNL)
INTEGER*4 IT,NPROJ
INTEGER*4 IB,IBL
INTEGER*4 IIRP,IIPL,ICDTT,ICDTTR
INTEGER*4 ICSUM,ICSUMS,ICDI,ISPA,ISPAT
INTEGER*4 KK0,KK0R,KK2,KKR,KDEE1,KDEE2
INTEGER*4 NE,NA,NG,NA1,NG1
INTEGER*4 JMAX,LJ,INEL,IH,IH1,IY,IY2,IY3
INTEGER*4 JL,KK1,IVMIN,IVMAX,NPA,IREC1,IREC,MAXA,NALL,NSA,KIS
INTEGER*4 IA,IAG,IAGS,IG,IESP,IESLOG
INTEGER*4 IPOT,IPOTR,IRL,ICDIR,ICSBR,ICSUMR,KOI,IGG,KIST
INTEGER*4 JRT,IESPT,IP,I1,IPB
INTEGER*4 IPB1,KIB,IPT,IE,IERLOG,IAGB,KIT,IMA,IIM
INTEGER*4 im1,im2,im3,IG2,ies,ias
INTEGER*4 JE,JA,JG,JTJ,JTK,JTL
INTEGER*4 NLayers
C REAL Variables
REAL*8 CVMGT
REAL*8 X(64),Y(64),Z(64),E(64),PL(64) ,COSX(64),COSY(64),COSZ(64)
& ,SINE(64)
REAL*8 EPS(64),DEN(64),DEE(64),DENS(64),DEES(64) ,CX(64),CY(64)
& ,CZ(64),SX(64),X1(64),ASIGT(64),EM(64)
REAL*8 EX1(64),EX2(64),EX3(64),P(64),TAU(64),EX4(64) ,B(64),R(64)
& ,C2(64),S2(64),CT(64),ST(64),V(64),V1(64) ,CPHI(64),SPHI(64)
& ,CPSI(64),SPSI(64),TAUPSI(64) ,ENUCL(64),EINEL(64),ENUCL2(64)
& ,EINEL2(64)
REAL*8 ER(2000,2),XR(2000,2),YR(2000,2),ZR(2000,2) ,CSXR(2000,2)
& ,CSYR(2000,2),CSZR(2000,2),TAUR(2000) ,SNXR(2000,2)
& ,CPSIR(2000,2),SPSIR(2000,2),CPHIR(2000,2) ,SPHIR(2000,2)
& ,TAUPSR(2000,2)
REAL*8 EPSR(2000),T(2000),TS(2000),DEER(2000),DEERS(2000)
& ,PR(2000),BR(2000),EX1R(2000),EX2R(2000),EX3R(2000),CTR(2000)
& ,STR(2000),ASIGTR(2000),EX4R(2000) ,X2(2000),RR(2000)
& ,VR(2000) ,V1R(2000),CXR(2000),CYR(2000),CZR(2000) ,SXR(2000)
& ,C2R(2000),S2R(2000),CUR(2000)
REAL*8 RIRP(0:MAXD1),CASMOT(MAXD),PHON(MAXD),DENT(MAXD),ION(MAXD)
& ,DMGN(MAXD),CASMOTR(MAXD),PHONR(MAXD),DENTR(MAXD),IONR(MAXD)
& ,DMGNR(MAXD),ELGD(MAXD),ELGDR(MAXD)
REAL*8 ELE(MAXD,MAXNL5),ELI(MAXD,MAXNL5),ELP(MAXD,MAXNL5)
& ,ELD(MAXD,MAXNL5) ,ELET(MAXNL5),ELIT(MAXNL5),ELPT(MAXNL5)
& ,ELDT(MAXNL5),ELER(MAXD,MAXNL5),ELIR(MAXD,MAXNL5),ELPR(MAXD
& ,MAXNL5),ELDR(MAXD,MAXNL5) ,ELETR(MAXNL5),ELITR(MAXNL5)
& ,ELPTR(MAXNL5),ELDTR(MAXNL5)
REAL*8 AI(20),RKADB(20),RKADT(20) ,RKADS(20),RKADST(20)
& ,RKADSJ(20,MAXNLm15),RKADSL(20,MAXNL),RKDSTJ(20,MAXNLm15)
& ,RKDSTL(20,MAXNL)
REAL*8 EBSP(MAXNL5),ESPAL(MAXNL) ,SPY(MAXNL5),SPE(MAXNL5)
& ,REY(MAXNL5),EMSP(MAXNL5),ESPIP(MAXNL5),ESPIS(MAXNL5)
& ,ESPOP(MAXNL5),ESPOS(MAXNL5) ,RIP(MAXNL5),RIS(MAXNL5)
& ,ROP(MAXNL5),ROS(MAXNL5) ,REIP(MAXNL5),REIS(MAXNL5)
& ,REOP(MAXNL5),REOS(MAXNL5),ESPMIP(MAXNL5),ESPMIS(MAXNL5)
& ,ESPMOP(MAXNL5),ESPMOS(MAXNL5) ,RIPJ(MAXNL5),RISJ(MAXNL5)
& ,ROPJ(MAXNL5),ROSJ(MAXNL5) ,REIPJ(MAXNL5),REISJ(MAXNL5)
& ,REOPJ(MAXNL5),REOSJ(MAXNL5)
REAL*8 ETSP(MAXNL5),ESPALT(MAXNL) ,SPYT(MAXNL5),SPET(MAXNL5)
& ,REYT(MAXNL5),EMSPT(MAXNL5),ESPIPT(MAXNL5),ESPIST(MAXNL5)
& ,ESPOPT(MAXNL5),ESPOST(MAXNL5),RIPT(MAXNL5),RIST(MAXNL5)
& ,ROPT(MAXNL5),ROST(MAXNL5) ,REIPT(MAXNL5),REIST(MAXNL5)
& ,REOPT(MAXNL5),REOST(MAXNL5) ,ESPMIPT(MAXNL5),ESPMIST(MAXNL5)
& ,ESPMOPT(MAXNL5),ESPMOST(MAXNL5)
REAL*8 SPEM(MAXNL5),SPE2S(MAXNL5),SPE3S(MAXNL5),SPE4S(MAXNL5)
& ,SPE5S(MAXNL5) ,SPE6S(MAXNL5),VSPE(MAXNL5),SSPE(MAXNL5)
& ,GSPE(MAXNL5),BSPE(MAXNL5)
REAL*8 SPE1SL(MAXNL5),SPE2SL(MAXNL5),SPE3SL(MAXNL5),SPE4SL(MAXNL5)
& ,SPE5SL(MAXNL5),SPE6SL(MAXNL5)
REAL*8 ELOG(75),EMA(62,22),EABL(75)
REAL*8 EMAT(62,22),EATL(75),EASL(75,MAXNLm15),EASTL(75,MAXNLm15)
REAL*8 FG(128),FFG(64)
REAL*8 XX(MAXNL),DX(MAXNL),RHO(MAXNL)
REAL*8 Z2(MAXNL),M2(MAXNL),LM(MAXNL)
REAL*8 PDMAX(MAXNL),ARHO(MAXNL),AM(MAXNL)
REAL*8 FM(MAXNL),EPS0(MAXNL),ASIG(MAXNL)
REAL*8 K2(MAXNL),CK(MAXNL),KLM1(MAXNL)
REAL*8 SB(MAXNL),DLI(MAXNL)
REAL*8 UpTiefe,LowTiefe
C ZT - atomin numbers, MT - mass numbers (amu), CO - concentration (stoichiometry)
C SBE - surface binding energy, ED - displacement energy, BE - bulk binding energy
C COM - ??
REAL*8 ZT(MAXNL,MAXEL),MT(MAXNL,MAXEL),CO(MAXNL,MAXEL)
& ,SBE(MAXNL,MAXEL),ED(MAXNL,MAXEL),BE(MAXNL,MAXEL),
& COM(MAXEL,MAXNL)
REAL*8 MU(MAXNL5,MAXNL5),EC(MAXNL5,MAXNL5),A(MAXNL5,MAXNL5)
& ,F(MAXNL5,MAXNL5) ,KL(MAXNL5,MAXNL5),KOR(MAXNL5,MAXNL5)
& ,KLM(MAXNL,MAXNL5)
REAL*8 MU1(MAXNL5),EC1(MAXNL5),A1(MAXNL5),F1(MAXNL5),KL1(MAXNL5)
& ,KOR1(MAXNL5) ,DI(MAXNL5),EP(MAXNL5),ZZ(MAXNL5),TM(MAXNL5)
C CH1,2,3,4,5 are values of A-1,2,3,4,5 of the Ziegler tables
REAL*8 CH1(MAXNL,MAXEL),CH2(MAXNL,MAXEL),CH3(MAXNL,MAXEL)
& ,CH4(MAXNL,MAXEL),CH5(MAXNL,MAXEL)
REAL*8 CHM1(MAXNL)
REAL*8 SM(64),SH(64,MAXEL)
REAL*8 FIESB(MAXNLm15),SEESB(MAXNLm15),THESB(MAXNLm15)
& ,FOESB(MAXNLm15) ,SGMESB(MAXNLm15),DFIESB(MAXNLm15)
& ,DSEESB(MAXNLm15),DTHESB(MAXNLm15)
REAL*8 pi,c,E0,de,alfa,z1,rtheta,cu,enot,esb,est,esp
REAL*8 E2,AB,FP,AN
REAL*8 Esig,Epar
REAL*8 E0keV,EsigkeV
c for part. reflec. coeff. calculation
REAL*8 epsilon, prcoeff,PRC(6)
REAL*8 cossin,rphi,vanlb,vailb,vanlt,vailt,phip,ta,ta2
REAL*8 exir,phipr,u,rq,es,vanli,vaili
REAL*8 M1,MOT
REAL*8 ET,PLST,PL2ST,PL3ST
REAL*8 PL4ST,PL5ST,PL6ST
REAL*8 SEM,TDMENR
REAL*8 TION,TIONR,TDENT,TDENTR
REAL*8 TELGD,TPHON,TCASMO,TELGDR
REAL*8 TPHONR,TRIRP,TDMGN,TDMGNR
REAL*8 ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM
REAL*8 EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM,EB
REAL*8 EB1SUL,EB2SUL,EB3SUL,EB4SUL,EB5SUL,EB6SUL
REAL*8 EINELB,EIL2B
REAL*8 EXI,exi1s,exi2s,exi3s,exi4s,exi5s,exi6s,exiq,exic
REAL*8 cossq,cosst,coss1s,coss2s,coss3s,coss4s,coss5s,coss6s
REAL*8 X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,XSUM
REAL*8 R2SUM,R3SUM,R4SUM,R5SUM,R6SUM,RSUM
REAL*8 PL2SUM,PL3SUM,PL4SUM,PL5SUM,PL6SUM,PLSUM
REAL*8 ENL2B,ENUCLB,EINELI,EIL2I
REAL*8 ENUCLI,ENL2I
REAL*8 PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB
REAL*8 EELWC,EELWC2,EELWC3,EELWC4,EELWC5,EELWC6
REAL*8 EIL,EIL2,EIL3,EIL4,EIL5,EIL6
REAL*8 EPL,EPL2,EPL3,EPL4,EPL5,EPL6
REAL*8 EEL,EEL2,EEL3,EEL4,EEL5,EEL6
REAL*8 EN2LT
REAL*8 EMX,ESPAT,ESPA
REAL*8 ALPHA,ALPHASIG,ALPHAPAR
REAL*8 EF,SHEATH,ERC,RI,RI2,RI3,X0,RD,CW,CA
REAL*8 DA,DG,DGI,BW,DAW,DGW,E0DE,DGIK,PI2,ABC
REAL*8 TT,HLM,HLMT
REAL*8 SU1,SU2,SUR,SU,SUT1,SUT2,SUTR,SUT
REAL*8 CALFA,SALFA,SFE,XC,RT,TI,SINA
REAL*8 ZARG,VELC,RA,RR1,FR,FR1,Q,FE(64)
REAL*8 ROCINV,SQE,CC,AA,FF,DELTA,DEEOR,DEEORR,DELR,FHE,DEL
REAL*8 G,DEN2,DEN3,DEE2,DEE3,DEWC,DEWC2,DEWC3
REAL*8 TAR,TAR2,RRR1,T1,TEL,TR,TR1,EI,ENOR,ACS,AC
REAL*8 SPE2,SPE3,SPE2L,SPE3L
REAL*8 ENORT,ESPT,EXIRT,EINELT,EIL2T
REAL*8 PL2,PL3,XQ,XQ3,RQW,RQ3,ENO,ESQ,ES3
REAL*8 ESQL,ES3L,PLQB,PL3B,PLQT,PL3T
REAL*8 ETQ,ET3,ENUCLT,ENL2T
REAL*8 RA1,ALPHAM,EMV,EIM
REAL*8 CSUM,CSUMS,CSUMR,AVCSUM,AVCSMS,AVCDIS
REAL*8 AVNLI,SIGNLI,DFINLI,SIGILI,DFIILI,AVILI
REAL*8 TIT,TE,TMEANR,EMEANT,AVNLT,SIGNLT,DFINLT,AVILT
REAL*8 TN,SIGILT,DFIILT
REAL*8 FIX0,SEX,THX,FOX,FIX,SIX,SIGMAX,DFIX0,DSEX,DTHX
REAL*8 FIR0,SER,THR,FOR,FIR,SIR,SIGMAR,DFIR0,DSER,DTHR
REAL*8 FIP0,SEP,THP,FOP,FIP,SIP,SIGMAP,DFIP0,DSEP,DTHP
REAL*8 FIE0,SEE,THE,FOE,FIE,SIE,SIGMAE,DFIE0,DSEE,DTHE
REAL*8 FIW0,SEW,THW,FOW,FIW,SIW,SIGMAW,DFIW0,DSEW,DTHW
REAL*8 FII0,SEI,THI,FOI,FII,SII,SIGMAI,DFII0,DSEI,DTHI
REAL*8 FIS0,SES,THS,FOS,FIS,SIS,SIGMAS,DFIS0,DSES,DTHS
REAL*8 FIB0,SEB,THB,FOB,FIB,SIB,SIGMAB,DFIB0,DSEB,DTHB
REAL*8 FIPB0,SEPB,THPB,FOPB,FIPB,SIPB,SIGMPB,DFIPB0,DSEPB,DTHPB
REAL*8 FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT
REAL*8 FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT,DFIPT0,DSEPT,DTHPT
REAL*8 FIES0,SEES,THES,FOES,FIES,SIES,SIGMES,DFIES0,DSEES,DTHES
REAL*8 FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,SIGMSL,DFIESL,DSEESL
& ,DTHESL
REAL*8 X1SD,X2SD,X3SD,X4SD,X5SD,X6SD
REAL*8 ACSUMR,ACDISR,ACSBER,ACSUR,ACDIR,ACSBR
REAL*8 ACDR11,ACDR12,ACDR21,ACDR22
REAL*8 D1,D2,Dmid,RN,RE,EMEANR,EMEAN,AVEB,AVNLB,SIGNLB
REAL*8 DFINLB,AVILB,SIGILB,DFIILB,TEMP,TEMPNH
REAL*8 EB1B,EB2B,EB3B,EB4B,EB5B,EB6B
REAL*8 EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL
REAL*8 EBSP1,EBSP2,EBSP3,EBSP4,EBSP5,EBSP6
REAL*8 EBSP1L,EBSP2L,EBSP3L,EBSP4L,EBSP5L,EBSP6L
REAL*8 PL1S,PL2S,PL3S,PL4S,PL5S,PL6S
REAL*8 YH,HN,CST,BI,BIL,YSP,YSPL,EEE
REAL*4 random(1),ran2(2)
C CHARACTER Variables
CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2
CHARACTER*60 filein,fileout
CHARACTER*64 innam,outnam,rgenam,errnam,seqnam
CHARACTER*4 inext,outext,rgeext,errext,seqext
CHARACTER errcom*72
CHARACTER month_start*4,month_stop*4,day_start*2,day_stop*2
CHARACTER year_start*4,year_stop*4,hour_start*2,hour_stop*2
CHARACTER min_start*2,min_stop*2,sec_start*2,sec_stop*2
CHARACTER*50 chem(MAXNL)
COMMON /A/ M1,VELC,ZARG
COMMON /B/ TI,SHEATH,CALFA
C=======================================================================
C STATIC INITIALISATION
C=======================================================================
DATA PI/3.14159265358979D0/, ICW/100/, E2/14.399651D0/
DATA AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/
DATA inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/
DATA errext/'.err'/,seqext/'.seq'/
DATA ET/0.D0/,PLST/0.D0/,PL2ST/0.D0/,PL3ST/0.D0/
DATA PL4ST/0.D0/,PL5ST/0.D0/,PL6ST/0.D0/
DATA KIT/0/,KIS/0/,KIST/0/,NALL/0/,NPA/0/,NSA/0/,KIB/0/,MAXA/0/
DATA SEM/0.D0/,IT/0/,NPROJ/0/,NREC1/0/,NREC2/0/
DATA NH/0/,IB/0/,IBL/0/,NJ/MAXNL*0/,NLI/0,MAXNL*0/
DATA DLI/MAXNL*0.D0/
DATA tryE/0/,negE/0/
DATA IIRP/0/,IIPL/0/,ICDTT/0/,ICDTTR/0/,TDMENR/0.D0/
DATA TION/0.D0/,TIONR/0.D0/,TDENT/0.D0/,TDENTR/0.D0/
DATA TELGD/0.D0/,TPHON/0.D0/,TCASMO/0.D0/,TELGDR/0.D0/
DATA TPHONR/0.D0/,TDMENR/0.D0/,TRIRP/0.D0/,TDMGN/0.D0/
DATA TDMGNR/0.D0/
DATA ET2SUM/0.D0/,ET3SUM/0.D0/,ET4SUM/0.D0/,ET5SUM/0.D0/
DATA ET6SUM/0.D0/
DATA EB2SUM/0.D0/,EB3SUM/0.D0/,EB4SUM/0.D0/,EB5SUM/0.D0/
DATA EB6SUM/0.D0/,EB/0.D0/
DATA EB1SUL/0.D0/,EB2SUL/0.D0/,EB2SUL/0.D0/,EB3SUL/0.D0/
DATA EB4SUL/0.D0/,EB5SUL/0.D0/,EB6SUL/0.D0/
DATA EINELB/0.D0/,EIL2B/0.D0/
DATA exi1s/0.D0/,exi2s/0.D0/,exi3s/0.D0/,exi4s/0.D0/
DATA exi5s/0.D0/,exi6s/0.D0/
DATA KADB/20*0/
DATA coss1s/0.D0/,coss2s/0.D0/,coss3s/0.D0/,coss4s/0.D0/
DATA coss5s/0.D0/,coss6s/0.D0/
DATA MEAB/MAXD2meab*0/,MEABL/1575*0/,MAGB/1364*0/
DATA MEPB/MAXD2p2*0/,MEAGB/MAXD2meagb*0/,EMA/1364*0.D0/
DATA X2SUM/0.D0/,X3SUM/0.D0/,X4SUM/0.D0/,X5SUM/0.D0/
DATA X6SUM/0.D0/,XSUM/0.D0/
DATA R2SUM/0.D0/,R3SUM/0.D0/,R4SUM/0.D0/,R5SUM/0.D0/
DATA R6SUM/0.D0/,RSUM/0.D0/
DATA PL2SUM/0.D0/,PL3SUM/0.D0/,PL4SUM/0.D0/,PL5SUM/0.D0/
DATA PL6SUM/0.D0/,PLSUM/0.D0/
DATA ENL2B/0.D0/,ENUCLB/0.D0/,EINELI/0.D0/,EIL2I/0.D0/
DATA ENUCLI/0.D0/,ENL2I/0.D0/
C BUGFIX: transmitted projectile loss accumulators were read before
C guaranteed assignment when at least one projectile transmits.
DATA ENUCLT/0.D0/,ENL2T/0.D0/,EINELT/0.D0/,EIL2T/0.D0/
DATA PLSB/0.D0/,PL2SB/0.D0/,PL3SB/0.D0/,PL4SB/0.D0/
DATA PL5SB/0.D0/,PL6SB/0.D0/
DATA EELWC/0.D0/,EELWC2/0.D0/,EELWC3/0.D0/,EELWC4/0.D0/
DATA EELWC5/0.D0/,EELWC6/0.D0/
DATA EIL/0.D0/,EIL2/0.D0/,EIL3/0.D0/,EIL4/0.D0/
DATA EIL5/0.D0/,EIL6/0.D0/
DATA EPL/0.D0/,EPL2/0.D0/,EPL3/0.D0/,EPL4/0.D0/
DATA EPL5/0.D0/,EPL6/0.D0/
DATA EEL/0.D0/,EEL2/0.D0/,EEL3/0.D0/,EEL4/0.D0/
DATA EEL5/0.D0/,EEL6/0.D0/
DATA ENUCL/64*0.D0/,EN2LT/0.D0/,TAUPSI/64*0.D0/
DATA EINEL/64*0.D0/,CASMOT/MAXD*0.D0/,DENT/MAXD*0.D0/
DATA DMGN/MAXD*0.D0/,ION/MAXD*0.D0/,PHON/MAXD*0.D0/
DATA PHONR/MAXD*0.D0/
DATA ELGD/MAXD*0.D0/,ELGDR/MAXD*0.D0/
DATA ICDT/MAXD*0/,ICDTR/MAXD*0/
DATA ICDR/MAXDNL5*0/,ICDIRN/MAXDNL5*0/,IONR/MAXD*0.D0/
DATA DENTR/MAXD*0.D0/,DMGNR/MAXD*0.D0/
DATA IPL/MAXD*0/,IPLB/MAXD*0/,IPLT/MAXD*0/,IRP/MAXD2*0/
DATA IRPL/MAXNL*0/
DATA ICDJT/MAXNL5*0/,ICDJTR/MAXNL5*0/,ICDITR/MAXNL5*0/
DATA ICD/MAXDNL5*0/,ELP/MAXDNL5*0.D0/,ELD/MAXDNL5*0.D0/
DATA ELE/MAXDNL5*0.D0/,ELI/MAXDNL5*0.D0/
DATA ICDIRI/MAXNL5p2*0/
DATA ICSUM/0/,ICSUMS/0/,ICDI/0/,ISPA/0/,ISPAT/0/
C BUGFIX: recoil counters are used in IRL output and incremented
C later; they must not start from compiler-dependent garbage.
DATA ICDIR/0/,ICSBR/0/,ICSUMR/0/
DATA Z2/MAXNL*0.D0/,M2/MAXNL*0.D0/
DATA KLM1/MAXNL*0.D0/,CHM1/MAXNL*0.D0/
DATA SB/MAXNL*0.D0/,KLM/MAXNLp25*0.D0/
DATA ME/5000*0/,EMX/0.D0/,ESPAT/0.D0/,ESPA/0.D0/
DATA IBSP/MAXNL5*0/,IBSPL/MAXNL5*0/,EBSP/MAXNL5*0.D0/
DATA ISPAL/MAXNL*0/
DATA ITSP/MAXNL5*0/,ETSP/MAXNL5*0.D0/
DATA ESPAL/MAXNL*0.D0/,ESPALT/MAXNL*0.D0/
DATA ISPALT/MAXNL*0/
DATA SPE2S/MAXNL5*0.D0/,SPE3S/MAXNL5*0.D0/,SPE4S/MAXNL5*0.D0/
DATA SPE5S/MAXNL5*0.D0/,SPE6S/MAXNL5*0.D0/
DATA SPE1SL/MAXNL5*0.D0/,SPE2SL/MAXNL5*0.D0/,SPE3SL/MAXNL5*0.D0/
DATA SPE4SL/MAXNL5*0.D0/,SPE5SL/MAXNL5*0.D0/,SPE6SL/MAXNL5*0.D0/
DATA ELET/MAXNL5*0.D0/,ELPT/MAXNL5*0.D0/,ELDT/MAXNL5*0.D0/
DATA ELIT/MAXNL5*0.D0/
DATA ELETR/MAXNL5*0.D0/,ELITR/MAXNL5*0.D0/,ELPTR/MAXNL5*0.D0/
DATA ELDTR/MAXNL5*0.D0/
DATA Epar/0.D0/
C part. refl. coeff. from Thomas et al.
DATA PRC/0.825D0,21.41D0,8.606D0,0.6425D0,1.907D0,1.927D0/
DATA number_in_layer/MAXNL*0/
C=======================================================================
C Variable dictionary (comment-only, fixed-form Fortran)
C
C Indices and counters
C
C I layer index unless used as a local loop counter
C J element index inside layer I
C K global element/species index over all layers
C IV active primary-projectile slot, 1..64
C IREC recoil-particle slot, 1..2000
C IP/IPB/IPT depth-bin indices for general/back/transmission use
C IE,IA,IG energy, polar-angle and azimuth-angle histogram indices
C
C Main dimensions
C
C MAXD number of depth bins
C MAXNL maximum number of target layers
C MAXEL maximum number of elements in one layer
C MAXNL5 maximum total layer-specific species, MAXNL*MAXEL
C
C Main run controls
C
C NH number of incident projectiles to simulate
C E0 initial projectile energy
C ESIG initial energy spread
C ALPHA incident angle
C X0 initial depth/start position
C CW depth-bin width
C L number of target layers used in this run
C
C Projectile
C
C Z1 projectile atomic number
C M1 projectile mass
C E(IV) current projectile energy
C X,Y,Z current projectile position components
C PL(IV) current projectile path length
C COSX,COSY,COSZ current projectile direction cosines
C
C Recoils
C
C ER recoil energy
C XR,YR,ZR recoil position components
C CSXR,CSYR,CSZR recoil direction cosines
C JJR recoil species index information
C LRR recoil layer index information
C INOUT recoil status / surface-crossing bookkeeping
C
C Target geometry and composition
C
C CHEM(I) layer name or chemical formula read from input
C DX(I) thickness of layer I
C RHO(I) mass density of layer I
C NJ(I) number of elements in layer I
C JT(I) cumulative species index up to layer I
C ZT(I,J) atomic number of element J in layer I
C MT(I,J) mass of element J in layer I
C CO(I,J) concentration fraction of element J in layer I
C Z2,M2 concentration-weighted layer atomic number and mass
C
C Binding, displacement and stopping input
C
C SBE(I,J) surface binding energy
C ED(I,J) displacement energy
C BE(I,J) bulk binding energy
C CH1-CH5 Ziegler stopping coefficients for each element
C CHM1(I) layer-level stopping coefficient helper
C
C Layer-specific species tables
C
C ZZ(K) atomic number for global species K
C TM(K) mass for global species K
C DI(K) density/concentration helper for species K
C EP(K) energy parameter for species K
C MU,EC,A,F,KL,KOR pair tables for projectile/target species
C MU1,EC1,A1,F1,KL1,KOR1 pair tables for incident projectile
C
C Energy loss and damage profiles
C
C ENUCL,EINEL projectile nuclear and electronic losses
C ENUCL2,EINEL2 squared-loss accumulators
C CASMOT,DENT,ION,PHON,DMGN depth-resolved projectile damage
C CASMOTR,DENTR,IONR,PHONR,DMGNR recoil-resolved damage
C ELE,ELI,ELP,ELD projectile energy-loss profiles by species
C ELER,ELIR,ELPR,ELDR recoil energy-loss profiles by species
C
C Range, reflection and transmission statistics
C
C IRP,RIRP implanted-projectile range histograms
C IPL,IPLB,IPLT path-length histograms
C IB,IT numbers of backscattered and transmitted projectiles
C ET sum of transmitted projectile energies
C EB sum of backscattered projectile energies
C
C Sputtering and recoil statistics
C
C IBSP,ITSP sputtered species counters, backward/transmission
C EBSP,ETSP sputtered species energy sums
C SPY,SPYT sputter yields per species
C REY,REYT recoil emission yields per species
C KADS/KADST sputtered angular counters
C KDSTJ/KDSTL transmitted sputtering angular counters
C
C Histogram arrays used in output
C
C MEAB,MEAT energy spectra for backscattered/transmitted ions
C MAGB,MAGT angular spectra for backscattered/transmitted ions
C MEAS,MEAST energy spectra for sputtered species
C MAGS,MAGST angular spectra for sputtered species
C EMA,EMAT accumulated energies in angular-energy bins
C
C File names and extensions
C
C FILEIN run basename supplied on the command line
C INNAM input file name, FILEIN//'.inp'
C OUTNAM main output file name, FILEIN//'.out'
C RGENAM range output file name, FILEIN//'.rge'
C ERRNAM error/log file name, FILEIN//'.err'
C=======================================================================
C=======================================================================
C COMMAND LINE AND FILE NAMES
C=======================================================================
CALL getarg(1, filein)
C Require an explicit run basename for input/output naming.
if (filein.eq.'') then
WRITE(*,*) 'Usage: trimspNL <run_basename>'
STOP 1
endif
fileout = filein
innam=trim(filein)//inext
outnam=trim(fileout)//outext
rgenam=trim(fileout)//rgeext
errnam=trim(fileout)//errext
seqnam=trim(fileout)//seqext
write (*,*) innam
C LMAX is maximum number of layers and JMAX is maximum number of
C elements per layer.
JMAX=MAXEL
C=======================================================================
C READ GUI-GENERATED INPUT FILE
C=======================================================================
C This part reads the input file (new format).
C The JavaScript/Electron frontend writes this sequential block layout,
C so the READ order below must stay in sync with CreateInpFile().
C Input reading has been moved to READ_GUI_INPUT so the
C GUI file layout is isolated from the transport code.
CALL READ_GUI_INPUT(INNAM,ERRNAM,Z1,M1,E0,ESIG,ALPHA,ALPHASIG,EF,
& ESB,SHEATH,ERC,NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,
& KDEE2,IPOT,IPOTR,IRL,NLayers,CHEM,DX,RHO,CK,NJ,ZT,MT,CO,SBE,
& ED,BE,CH1,CH2,CH3,CH4,CH5)
C=======================================================================
C OPEN OUTPUT FILES AND START CLOCK
C=======================================================================
C open statement for output files, removed from line 2449 ff to here
OPEN(UNIT=21,FILE=outnam)
6001 OPEN(UNIT=22,FILE=rgenam,STATUS='replace')
WRITE(21,1000) TRIMSP_VERSION
1000 FORMAT(1H1/,6X,'* TrimSPNL v',A,' *')
C Get simulation start time
CALL TimeStamp(day_start,month_start,year_start,hour_start
& ,min_start,sec_start,seconds_start_total)
WRITE(21,*)
WRITE(21,1050)day_start,month_start,year_start,hour_start
& ,min_start,sec_start
WRITE(*,1001) TRIMSP_VERSION
1001 FORMAT(1X,'Version ',A)
WRITE(*,1050)day_start,month_start,year_start,hour_start
& ,min_start,sec_start
1050 FORMAT(1x,' Start: ',A2,'.',A4,1x,A4,1x,A2
& ,':',A2,':',A2)
C=======================================================================
C OUTPUT BINNING CONSTANTS
C=======================================================================
C SET INTERVAL CONSTANTS FOR OUTPUT
DE = 1.D0
DA = 3.D0
DG = 3.D0
DGI = 15.D0
NE = IDINT(100.D0/DE + 2.00001D0)
NA = IDINT(90.D0/DA + 2.00001D0)
NG = IDINT(180.D0/DG + 2.00001D0)
NG = NA +NA -2
NGIK = IDINT(180.D0/DGI+ 0.001D0)
NE1 = NE -1
NA1 = NA -1
NG1 = NG -1
IF(E0.LT.0.) then
E0DE = 10.0D0/(DABS(E0)*DE)
else
E0DE = 100.0D0/(E0*DE)
endif
BW = 180.D0/PI
DAW = BW/DA
DGW = BW/DG
DGIK = BW/DGI
C=======================================================================
C TARGET SETUP AND PRECOMPUTED TABLES
C=======================================================================
C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS
PI2=2.D0*PI
ABC=AB*FP
L=NLayers
C Checks wether depth interval is an integer denominator of layer thickness or not.
C If not, the calculated number of implanted particles in each layer from bin width
C is wrong. Calculation form layer thickness should be correct.
C The rge files (stopping profiles) and summary files should be fine.
C CW is depth increment and DX(K) is the thickness of layer K
depth_interval_flag = 1
C DO K=1,L-1
C IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN
C If the thickness of layer K is not an integer multiples of depth increment
C depth_interval_flag = 0
C GO TO 44
C ENDIF
C ENDDO
C 44 CONTINUE
JT(1) = 0
JT(2) = NJ(1)
LJ= NJ(1)+NJ(2)
XX(1)=DX(1)
NJJ=0
C Loop over all defined layers
DO I=1,L
C BUGFIX: this loop iterates with I; using K here read an
C uninitialised/stale index and could test the wrong layer.
if(.not.EQUAL(DX(I)/CW-DBLE(IDINT(DX(I)/CW)),0.D0)) then
depth_interval_flag = 0
endif
if (I.ge.2) XX(I)=XX(I-1)+DX(I)
if (I.ge.3) then
JT(I)=JT(I-1)+NJ(I-1)
LJ=LJ+NJ(I)
endif
DO J=1,NJ(I)
Z2(I)=Z2(I)+CO(I,J)*ZT(I,J)
M2(I)=M2(I)+CO(I,J)*MT(I,J)
ENDDO
ARHO(I) = RHO(I)*AN/M2(I)
LM(I) = ARHO(I)**(-1.D0/3.D0)
ASIG(I) = LM(I)*ARHO(I)
PDMAX(I) = LM(I)/DSQRT(PI)
K2(I) = .133743D0*Z2(I)**(2.D0/3.D0)/DSQRT(M2(I))
AM(I) = CA*ABC*(Z2(I)**(-1.D0/3.D0))
FM(I) = AM(I)*M2(I)/(Z1*Z2(I)*E2*(M1+M2(I)))
EPS0(I) = FM(I)*E0
C For each layer calculate the following
do J=1,NJ(I)
ZZ(NJJ+J) = ZT(I,J)
TM(NJJ+J) = MT(I,J)
DI(NJJ+J) = ED(I,J)
EP(NJJ+J) = BE(I,J)
enddo
NJJ=NJJ+NJ(I)
COM(1,I) = CO(I,1)
DO J=1,NJ(I)-1
COM(J+1,I) = COM(J,I)+CO(I,J+1)
ENDDO
ENDDO
DO J = 1,LJ
MU1(J) = M1/TM(J)
EC1(J) = 4.D0*MU1(J)/((1.D0+MU1(J))*(1.D0+MU1(J)))
C KR-C (IPOT=1), MOLIERE (IPOT=2), ZBL POTENTIAL (IPOT=3)
A1(J) = CVMGT(CA*ABC*(ZZ(J)**(-1.D0/3.D0)),CA*ABC/(Z1**0.23D0
& +ZZ(J)**0.23D0),IPOT.LT.POT_ZBL)
F1(J) = A1(J)*TM(J)/(Z1*ZZ(J)*E2*(M1+TM(J)))
KL1(J) = 1.212D0*Z1**(7.D0/6.D0)*ZZ(J)/ ((Z1**(2.D0/3.D0)+ZZ(J)
& **(2.D0/3.D0))**1.5D0*DSQRT(M1))
ENDDO
IF(IPOT.EQ.POT_KRC) THEN
C KR-C POTENTIAL (IPOT=1)
DO J=1,LJ
KOR1(J) = 0.0389205D0*KL1(J)/(PI*A1(J)*A1(J))
ENDDO
ELSEIF (IPOT.EQ.POT_MOLIERE) THEN
C MOLIERE POTENTIAL (IPOT=2)
DO J=1,LJ
KOR1(J) = 0.045D0*KL1(J)/(PI*A1(J)*A1(J))
ENDDO
ELSEIF (IPOT.EQ.POT_ZBL) THEN
C ZBL POTENTIAL
DO J=1,LJ
KOR1(J) = 0.0203253D0*KL1(J)/(PI*A1(J)*A1(J))
ENDDO
ENDIF
DO I = 1,LJ
DO J = 1,LJ
MU(I,J) = TM(I)/TM(J)
EC(I,J) = 4.D0*MU(I,J)/((1.D0+MU(I,J))*(1.D0+MU(I,J)))
C KR-C , MOLIERE , ZBL POTENTIAL
A(I,J)= CVMGT(CA*ABC/(DSQRT(ZZ(I))+DSQRT(ZZ(J)))**(2.D0/3.D0
& ),CA*ABC/(ZZ(I)**0.23D0+ZZ(J)**0.23D0),
& IPOTR.LT.POT_ZBL)
C ZBL POTENTIAL
F(I,J) = A(I,J)*TM(J)/(ZZ(I)*ZZ(J)*E2*(TM(I)+TM(J)))
KL(I,J) = 1.212D0*ZZ(I)**(7.D0/6.D0)*ZZ(J)/ ((ZZ(I)**(2.D0
& /3.D0)+ZZ(J)**(2.D0/3.D0))**1.5D0*DSQRT(TM(I)))
ENDDO
ENDDO
IF (IPOTR.EQ.POT_KRC) THEN
C KR-C POTENTIAL (IPOTR=1)
DO I = 1,LJ
DO J = 1,LJ
KOR(I,J) = 0.0389205D0*KL(I,J)/(PI*A(I,J)*A(I,J))
ENDDO
ENDDO
ELSEIF (IPOTR.EQ.POT_MOLIERE) THEN
C MOLIERE POTENTIAL (IPOTR=2)
DO I = 1,LJ
DO J = 1,LJ
KOR(I,J) = 0.045D0*KL(I,J)/(PI*A(I,J)*A(I,J))
ENDDO
ENDDO
ELSEIF (IPOTR.EQ.POT_ZBL) THEN
C ZBL POTENTIAL (IPOTR=3)
DO I = 1,LJ
DO J = 1,LJ
KOR(I,J) = 0.0203253D0*KL(I,J)/(PI*A(I,J)*A(I,J))
ENDDO
ENDDO
ENDIF
DO LL=1,L
DO J=1,NJ(LL)
KLM1(LL) = KLM1(LL)+CO(LL,J)*KL1(J+JT(LL))*CK(LL)
CHM1(LL) = CHM1(LL)+CO(LL,J)*CH1(LL,J)
SB(LL) = SB(LL)+CO(LL,J)*SBE(LL,J)
ENDDO
ENDDO
DO I=1,LJ
DO LL = 1,L
DO J=1,NJ(LL)
KLM(LL,I) = KLM(LL,I)+CO(LL,J)*KL(I,J+JT(LL))
ENDDO
ENDDO
ENDDO
ALPHA = CVMGT( .001D0, ALPHA, EQUAL(ALPHA,0.D0))
ALPHA = CVMGT( 179.999D0, ALPHA, EQUAL(ALPHA,180.D0))
IF(ALPHA.GE.90.0.AND.X0.LE.0.0) GO TO 8881
GO TO 8882
8881 WRITE(6,8883)
8883 FORMAT(1X,'ERROR : IF ALPHA.GE.90. THEN IT MUST BE X0.LE.0.')
GO TO 8000
8882 CONTINUE
C SET CONSTANT DISTANCES
TT = XX(L)
INEL = 0
HLM = CVMGT( 0.D0, -.5D0*LM(1), INEL.EQ.0)
HLMT = CVMGT( TT, TT+.5D0*LM(L), INEL.EQ.0)
SU1 = PDMAX(1) + PDMAX(1)
SU2 = PDMAX(1)*(1.D0+KK0)
SUR = PDMAX(1)*(1.D0+KK0R)
SU = DMAX1(SUR,DMAX1(SU1,SU2))
SUT1 = TT + PDMAX(L) + PDMAX(L)
SUT2 = TT + PDMAX(L)*(1.D0+KK0)
SUTR = TT + PDMAX(L)*(1.D0+KK0R)
SUT = DMAX1(SUTR,DMAX1(SUT1,SUT2))
XC = CVMGT( X0, -SU, X0.GT.0.D0)
RT = TT-RD
IF(E0.GE.0.D0) GO TO 51
C
C=======================================================================
C INITIAL PROJECTILE BATCH
C=======================================================================
C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION
C
TI = -1.D0*E0
ZARG = DSQRT(TI/(M1*2.D0))
VELC = SHEATH/M1
C
C NUMBERS FOR VECTORIZED LOOPS
C
51 NUM = MIN( 64, NH)
IH = NUM
IH1 = NUM
C
C RANDOM START CONDITIONS
C
IY = IDINT(RI)
IY2 = IDINT(RI2)
IY3 = IDINT(RI3)
ISEED = IY
ISEED2 = IY2
ISEED3 = IY3
IF ( E0.GT.0.D0 ) GO TO 47
IF ( ALPHA.GE.0.D0 ) THEN
C MAXWELLIAN VELOCITY DISTRIBUTION
CALL VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,NUM)
DO IV=1,NUM
EMX = EMX+E(IV)
ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0))
C me(ne) is undefined here!
write (*,*) 'me=',me(ne)
me(ne) = me(ne)+1
ENDDO
GO TO 56
ELSE
CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM)
DO IV=1,NUM
EMX = EMX+E(IV)
ENDDO
GO TO 56
ENDIF
47 CONTINUE
IF(EQUAL(Esig,0.D0)) THEN
C FIXED PROJECTILE ENERGY
DO IV=1,NUM
E(IV) = E0
ENDDO
ELSE
C GAUSSIAN ENERGY DISTRIBUTION
DO IV=1,NUM
5200 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0)
tryE = tryE+1
IF(Epar.LE.0.0D0) THEN
negE = negE+1
GO TO 5200
ENDIF
E(IV)=Epar
ENDDO
ENDIF
SFE = DMIN1(SB(1),SB(L))
IF ( ALPHA.GE.0.D0 ) THEN
IF(EQUAL(ALPHASIG,0.D0))THEN
C
C FIXED PROJECTILE ANGLE
C
ALFA = ALPHA /BW
CALFA = DCOS(ALFA)
SALFA = DSIN(ALFA)
DO IV=1,NUM
COSX(IV) = CALFA
COSY(IV) = SALFA
COSZ(IV) = 0.D0
SINE(IV) = COSY(IV)
ENDDO
ELSE
C
C 1-D GAUSSIAN DISTRIBUTION OF ANGLE
C
DO IV=1,NUM
CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar,
& CALFA,SALFA,BW)
COSX(IV) = CALFA
COSY(IV) = SALFA
COSZ(IV) = 0.D0
SINE(IV) = COSY(IV)
ENDDO
ENDIF
ELSEIF (EQUAL(ALPHA,-2.D0)) THEN
C
C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL)
C
DO IV=1,NUM
call ranlux(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = DSQRT(RTHETA)
SINE(IV) = DSQRT(1.D0-RTHETA)
COSY(IV) = SINE(IV)*DCOS(RPHI)
COSZ(IV) = SINE(IV)*DSIN(RPHI)
ENDDO
ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN
C
C RANDOM DISTRIBUTION
C
DO IV=1,NUM
call ranlux(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = 1.D0 -2.D0*RTHETA
SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV))
COSY(IV) = SINE(IV) *DSIN(RPHI)
COSZ(IV) = SINE(IV) *DCOS(RPHI)
ENDDO
ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.LE.0.D0) THEN
DO IV=1,NUM
call ranlux(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = 1.D0 -RTHETA
SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV))
COSY(IV) = SINE(IV) *DSIN(RPHI)
COSZ(IV) = SINE(IV) *DCOS(RPHI)
ENDDO
ENDIF
56 IF ( X0.GT.0.D0 ) GO TO 59
C
C EXTERNAL START
C
DO IV=1,NUM
SINA = SINE(IV)
COSX(IV) = DSQRT( ( E(IV)*COSX(IV)*COSX(IV)+ESB)/(E(IV)+ESB))
SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV))
COSY(IV) = COSY(IV) *SINE(IV) /SINA
COSZ(IV) = COSZ(IV) *SINE(IV) /SINA
E(IV) = E(IV) + ESB
ENDDO
C
C LOCUS OF FIRST COLLISION
C
59 JL = ISRCHFGT(L,XX(1),1,X0)
DO IV=1,NUM
call ranlux(random, 1)
RA = CVMGT(DBLE(random(1)),1.D0,X0.LE.0.0D0)
X(IV) = XC + LM(JL) *RA *COSX(IV)
Y(IV) = LM(JL) *RA *COSY(IV)
Z(IV) = LM(JL) *RA *COSZ(IV)
PL(IV) = CVMGT(0.D0,LM(JL)*RA,X0.LE.0.0)
LLL(IV) = JL
ENDDO
C
C=======================================================================
C PRIMARY PROJECTILE TRANSPORT LOOP
C=======================================================================
C PROJECTILE LOOP
C Each pass transports one projectile history. The state vectors
C X/Y/Z, E, direction cosines, current layer index, etc. are updated
C in-place until the particle stops, backscatters or transmits.
C
1 CONTINUE
NPROJ=NPROJ+1
DO IV=1,IH1
CX(IV)=COSX(IV)
CY(IV)=COSY(IV)
CZ(IV)=COSZ(IV)
SX(IV)=SINE(IV)
DEES(IV)=0.D0
DENS(IV)=0.D0
DEN(IV)=0.D0
ENDDO
KK1=KK0
C
C=======================================================================
C PROJECTILE COLLISION CALCULATION
C=======================================================================
C COLLISION LOOP (INCLUDES WEAK SIMULTANEOUS COLL. FOR KK1.LT.4)
C KK controls the treatment of weak simultaneous collisions.
C For every active projectile we choose a collision partner, sample an
C impact parameter, solve the scattering geometry for the selected
C interatomic potential and then update the projectile state.
C
DO 245 KK=KK1,0,-1
C
C CHOICE OF COLLISION PARTNERS
C
DO IV=1,IH1
call ranlux(random, 1)
JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1,
& DBLE(random(1)))+JT(LLL(IV))
EPS(IV)=E(IV)*F1(JJJ(IV))
ENDDO
DO IV=1,IH1
C
C RANDOM AZIMUTHAL ANGLE AND IMPACT PARAMETER
C
call ranlux(ran2, 2)
PHIP=PI2*DBLE(ran2(1))
CPHI(IV)=DCOS(PHIP)
SPHI(IV)=DSIN(PHIP)
P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(ran2(2))+KK)
C
C POSITION OF TARGET ATOM
C
X1(IV)=X(IV)-P(IV)*CPHI(IV)*SX(IV)
P(IV)=CVMGT(1.D10,P(IV),X1(IV).LT.0.D0.OR.X1(IV).GT.TT)
B(IV)=P(IV)/A1(JJJ(IV))
ENDDO
CALL SCOPY(IH1,B,1,R,1)
IVMIN=1
IVMAX=IH1
C
C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.)
C
IF(IPOT.NE.POT_KRC) GO TO 4101
C KRYPTON-CARBON POTENTIAL
104 DO IV=IVMIN,IVMAX
IF(R(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 104 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
R(IV)=0.00001D0
ENDIF
EX1(IV)=DEXP(-.278544D0*R(IV))
EX2(IV)=DEXP(-.637174D0*R(IV))
EX3(IV)=DEXP(-1.919249D0*R(IV))
RR1=1.D0/R(IV)
V(IV)=(.190945D0*EX1(IV)+.473674D0*EX2(IV)+.335381D0*EX3(IV))
& *RR1
FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV)
V1(IV)=-(V(IV)+.05318658408D0*EX1(IV)+.301812757276D0*EX2(IV)+
& .643679648869D0*EX3(IV))*RR1
FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0
Q=FR/FR1
R(IV)=R(IV)-Q
TEST(IV)=DABS(Q/R(IV)).GT.0.001D0
ENDDO
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 106
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 106
GO TO 104
106 DO IV=1,IH1
ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV))
SQE=DSQRT(DABS(EPS(IV)))
CC=(.235809D0+SQE)/(.126000D0+SQE)
AA=2.D0*EPS(IV)*(1.D0+(1.0144D0/SQE))*B(IV)**CC
FF=(DSQRT(AA*AA+1.)-AA)*((69350.D0+EPS(IV))/(83550.D0+EPS(IV)))
DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0)
C2(IV)=DMIN1(1.0D0,C*C)
S2(IV)=1.D0-(1.D0*C2(IV))
ENDDO
GO TO 4103
4101 IF(IPOT.NE.POT_MOLIERE) GO TO 4102
C MOLIERE POTENTIAL
C CALL MAGICMOL(C2(1),S2(1),B(1),R(1),EPS(1),IH1)
4104 DO IV=IVMIN,IVMAX
IF(R(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 4104 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
R(IV)=0.00001D0
ENDIF
EX1(IV)=DEXP(-.3D0*R(IV))
EX2(IV)=DEXP(-1.2D0*R(IV))
EX3(IV)=DEXP(-6.0D0*R(IV))
RR1=1.D0/R(IV)
V(IV)=(.35D0*EX1(IV)+.55D0*EX2(IV)+.10D0*EX3(IV))*RR1
FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV)
V1(IV)=-(V(IV)+.105D0*EX1(IV)+.66D0*EX2(IV)+.6D0*EX3(IV))*RR1
FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0
Q=FR/FR1
R(IV)=R(IV)-Q
TEST(IV)=DABS(Q/R(IV)).GT.0.001D0
ENDDO
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 4106
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 4106
GO TO 4104
4106 DO IV=1,IH1
ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV))
SQE=DSQRT(EPS(IV))
CC=(.009611D0+SQE)/(.005175D0+SQE)
AA=2.D0*EPS(IV)*(1.D0+(0.6743D0/SQE))*B(IV)**CC
FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPS(IV))/(10.D0+EPS(IV)))
DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0)
C2(IV)=DMIN1(1.0D0,C*C)
S2(IV)=1.D0-(1.D0*C2(IV))
ENDDO
GO TO 4103
4102 IF(IPOT.NE.POT_ZBL) GO TO 4103
C ZBL POTENTIAL
C CALL MAGICZBL(C2(1),S2(1),B(1),R(1),EPS(1),IH1)
5104 DO IV=IVMIN,IVMAX
IF(R(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 5104 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
R(IV)=0.00001D0
ENDIF
EX1(IV)=DEXP(-.20162D0*R(IV))
EX2(IV)=DEXP(-.4029D0*R(IV))
EX3(IV)=DEXP(-.94229D0*R(IV))
EX4(IV)=DEXP(-3.1998D0*R(IV))
RR1=1.D0/R(IV)
V(IV)=(.02817D0*EX1(IV)+.28022D0*EX2(IV)+.50986D0*EX3(IV)+
& .18175D0*EX4(IV))*RR1
FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV)
V1(IV)=-(V(IV)+.0056796354D0*EX1(IV)+.112900638D0*EX2(IV)+
& .4804359794D0*EX3(IV)+.58156365D0*EX4(IV))*RR1
FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0
Q=FR/FR1
R(IV)=R(IV)-Q
TEST(IV)=DABS(Q/R(IV)).GT.0.001D0
ENDDO
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 5106
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 5106
GO TO 5104
5106 DO IV=1,IH1
IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 5106 '
ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV))
SQE=DSQRT(EPS(IV))
CC=(.011615D0+SQE)/(.0071222D0+SQE)
AA=2.D0*EPS(IV)*(1.D0+(0.99229D0/SQE))*B(IV)**CC
FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPS(IV))/(14.813D0+EPS(IV)
& ))
DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0)
C2(IV)=DMIN1(1.0D0,C*C)
S2(IV)=1.D0-(1.D0*C2(IV))
ENDDO
4103 CONTINUE
C
C END OF MAGIC
C
DO IV=1,IH1
DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV)
IF(C2(IV).LT.1.D-10) THEN
WRITE(99,'(A50)')' C2 < 10^-10, C2,S2,DEN resettet '
C2(IV)=1.D-10
S2(IV)=1.D0-(1.D0*C2(IV))
DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV)
ENDIF
TAU(IV)=CVMGT(P(IV)*DSQRT(DABS(S2(IV)/C2(IV))),0.D0,KK.EQ.0)
TAU(IV)=DMIN1(TAU(IV),LM(LLL(IV)))
CT(IV)=C2(IV)+C2(IV)-1.D0
ST(IV)=DSQRT(DABS(1.D0-CT(IV)*CT(IV)))
CU=CT(IV)+MU1(JJJ(IV))
CU=CVMGT(CU,1.0D-8,DABS(CU).GE.1.0D-8)
TA=ST(IV)/CU
TA2=1.D0/DSQRT(DABS(1.D0+TA*TA))
CPSI(IV)=CVMGT(TA2,-TA2,CU.GT.0.D0)
SPSI(IV)=DABS(TA)*TA2
DEEOR=CVMGT(KOR1(JJJ(IV))*DSQRT(DABS(E(IV)))*EX1(IV),0.D0,
& KDEE1.EQ.STOP_OR.OR.KDEE1.EQ.STOP_MIXED)
DENS(IV)=DENS(IV)+DEN(IV)
DEES(IV)=DEES(IV)+DEEOR
ENDDO
C
C DETERMINATION OF NEW FLIGHT DIRECTIONS
C
CALL DIRCOS(COSX(1),COSY(1),COSZ(1),SINE(1),CPSI(1),SPSI(1)
& ,CPHI(1),SPHI(1),IH1)
245 CONTINUE
C
C END OF COLLISION LOOP
C
C=======================================================================
C PROJECTILE ENERGY LOSS AND DAMAGE ACCUMULATION
C=======================================================================
C INELASTIC ENERGY LOSS( 5 POSSIBILITIES)
C
DO IV=1,IH1
ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV))
TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV))
ENDDO
C KDEE1 uses STOP_LS, STOP_OR, STOP_MIXED, STOP_ICRU49,
C and STOP_ZIEGLER in this order.
GO TO(15,16,17,18,19),KDEE1
15 DO IV=1,IH1
DEE(IV)=CVMGT(0.D0,KLM1(LLL(IV))*ASIGT(IV)*DSQRT(E(IV)),X(IV)
& .LT.HLM.OR.X(IV).GT.HLMT)
ENDDO
GO TO 40
16 DO IV=1,IH1
DEE(IV)=DEES(IV)
ENDDO
GO TO 40
17 DO IV=1,IH1
DEE(IV)=CVMGT(DEES(IV),0.5D0*(KLM1(LLL(IV))*ASIGT(IV)
& * DSQRT(E(IV))+DEES(IV)),X(IV).LT.HLM.OR.X(IV).GT.HLMT)
ENDDO
GO TO 40
18 DO IV=1,IH1
SM(IV)=0.D0
EM(IV)=E(IV)*0.001D0/M1
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SH(IV,J)=CVMGT(CH1(LLL(IV),J)*DSQRT(EM(IV)) ,CH2(LLL(IV),J)
& *EM(IV)**0.45D0*(CH3(LLL(IV),J)/EM(IV)) *DLOG(DABS(1.D0
& +(CH4(LLL(IV),J)/ EM(IV))+CH5(LLL(IV),J)*EM(IV)))
& /(CH2(LLL(IV),J)*EM(IV)**0.45D0+(CH3(LLL(IV),J)/EM(IV))
& *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/ EM(IV))+CH5(LLL(IV),J)
& *EM(IV)))) ,EM(IV).LT.10.D0)
ENDDO
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J)
ENDDO
ENDDO
DO IV=1,IH1
DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10
& .D0)
ENDDO
DO IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)*CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR
& .X(IV).GT.HLMT)
ENDDO
GO TO 40
19 FHE=CVMGT(1.3333D0,1.D0,M1.LT.4.00D0)
DO IV=1,IH1
SM(IV)=0.D0
EM(IV)=E(IV)*0.001D0*FHE
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SH(IV,J)=CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)* (CH3(LLL(IV)
& ,J)/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J)
& /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0))
& /(CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)+ (CH3(LLL(IV),J
& )/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J)
& /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0)))
ENDDO
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J)
ENDDO
ENDDO
DO IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)* CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV
& ).GT.HLMT)
ENDDO
40 CONTINUE
C
DO IV=1,IH1
DEL=DMAX1(1.0D-20,DENS(IV)+DEE(IV))
DENS(IV)=CVMGT(E(IV)*DENS(IV)/DEL,DENS(IV),DEL.GT.E(IV))
DEE(IV)=CVMGT(E(IV)*DEE(IV)/DEL,DEE(IV),DEL.GT.E(IV))
ENDDO
C
C INCREMENT OF DAMAGE, CASCADE AND PHONON ENERGY
C
DO 70 IV=1,IH1
I=MAX0(MIN0(IDINT(X1(IV)/CW+1.D0),100),1)
DENT(I)=DENT(I)+DENS(IV)
DMGN(I)=DMGN(I)+DEN(IV)
ION(I)=ION(I)+DEE(IV)
ELE(I,JJJ(IV))=ELE(I,JJJ(IV))+DEN(IV)
ELI(I,JJJ(IV))=ELI(I,JJJ(IV))+DEE(IV)
IF(DEN(IV).LE.DI(JJJ(IV))) GO TO 28
EPS(IV)=F1(JJJ(IV))*DEN(IV)
G=EPS(IV)+.40244D0*EPS(IV)**.75D0+3.4008D0*EPS(IV)**.16667D0
MOT=DEN(IV)/(1.D0+K2(LLL(IV))*G)
CASMOT(I)=CASMOT(I)+MOT
ELGD(I)=ELGD(I)+DEN(IV)
ELD(I,JJJ(IV))=ELD(I,JJJ(IV))+DEN(IV)
ICD(I,JJJ(IV))=ICD(I,JJJ(IV))+1
GO TO 70
28 PHON(I)=PHON(I)+DEN(IV)
ELP(I,JJJ(IV))=ELP(I,JJJ(IV))+DEN(IV)
70 CONTINUE
DO IV=1,IH1
ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV))))
ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1)))
ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0))
DEN2=DEN(IV)*DEN(IV)
DEN3=DEN2*DEN(IV)
EEL=EEL+DEN(IV)
EEL2=EEL2+DEN2
EEL3=EEL3+DEN3
EEL4=EEL4+DEN2*DEN2
EEL5=EEL5+DEN3*DEN2
EEL6=EEL6+DEN3*DEN3
DEE2=DEE(IV)*DEE(IV)
DEE3=DEE2*DEE(IV)
EIL=EIL+DEE(IV)
EIL2=EIL2+DEE2
EIL3=EIL3+DEE3
EIL4=EIL4+DEE2*DEE2
EIL5=EIL5+DEE3*DEE2
EIL6=EIL6+DEE3*DEE3
EPL=EPL+CVMGT(DEN(IV),0.D0,DEN(IV).LT.DI(JJJ(IV)))
EPL2=EPL2+CVMGT(DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV)))
EPL3=EPL3+CVMGT(DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV)))
EPL4=EPL4+CVMGT(DEN2*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV)))
EPL5=EPL5+CVMGT(DEN3*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV)))
EPL6=EPL6+CVMGT(DEN3*DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV)))
ENUCL(IV)=ENUCL(IV)+DENS(IV)
EINEL(IV)=EINEL(IV)+DEE(IV)
ENDDO
IF(KK0.EQ.0) GO TO 89
DO IV=1,IH1
DEWC=DENS(IV)-DEN(IV)
DEWC2=DEWC*DEWC
DEWC3=DEWC2*DEWC
EELWC=EELWC+DEWC
EELWC2=EELWC2+DEWC2
EELWC3=EELWC3+DEWC3
EELWC4=EELWC4+DEWC2*DEWC2
EELWC5=EELWC5+DEWC3*DEWC2
EELWC6=EELWC6+DEWC3*DEWC3
ENDDO
89 CONTINUE
IF(IRL.EQ.IRL_OFF) GO TO 27
C
C=======================================================================
C RECOIL GENERATION AND TRANSPORT
C=======================================================================
C VECTORIZED RECOIL LOOP
C
C TARGET RECOIL ATOM SECTION
C
C PRIMARY KNOCK-ON ATOMS
C
DO 6 IV=1,IH1
IF(DEN(IV).LE.ERC) GO TO 6
IF(X1(IV).GT.RD.AND.X1(IV).LT.RT) GO TO 6
NREC1=NREC1+1
ER(NREC1,1)=DEN(IV)-EP(JJJ(IV))
XR(NREC1,1)=X1(IV)
YR(NREC1,1)=Y(IV)-P(IV)*(SPHI(IV)*CZ(IV)-CPHI(IV)*CY(IV)*CX(IV)
& )/SX(IV)
ZR(NREC1,1)=Z(IV)+P(IV)*(SPHI(IV)*CY(IV)+CPHI(IV)*CX(IV)*CZ(IV)
& )/SX(IV)
CSXR(NREC1,1)=CX(IV)
CSYR(NREC1,1)=CY(IV)
CSZR(NREC1,1)=CZ(IV)
SNXR(NREC1,1)=SX(IV)
CPHIR(NREC1,1)=-CPHI(IV)
SPHIR(NREC1,1)=-SPHI(IV)
CT(IV)=DMIN1(CT(IV),.99999999D0)
TAR=ST(IV)/(1.D0-CT(IV))
TAR2=1./DSQRT(1.D0+TAR*TAR)
CPSIR(NREC1,1)=TAR2
SPSIR(NREC1,1)=TAR*TAR2
TAUPSR(NREC1,1)=0.D0
JJR(NREC1,1)=JJJ(IV)
KO(NREC1,JJR(NREC1,1),1)=1
INOUT(NREC1,1)=SIGN(1.D0,CX(IV))
NPA=NPA+1
6 CONTINUE
C
IF(NREC1.LT.NUM) GO TO 27
C
C START PROCESSING THE TARGET RECOIL ATOMS
C
83 CONTINUE
C
CALL DIRCOS(CSXR(1,1),CSYR(1,1),CSZR(1,1),SNXR(1,1)
1 ,CPSIR(1,1),SPSIR(1,1),CPHIR(1,1),SPHIR(1,1),NREC1)
C
C MOVE TARGET RECOIL ATOMS TO LIST 2
DO IREC1=1,NREC1
IREC=IREC1+NREC2
ER(IREC,2)=ER(IREC1,1)
XR(IREC,2)=XR(IREC1,1)
YR(IREC,2)=YR(IREC1,1)
ZR(IREC,2)=ZR(IREC1,1)
CSXR(IREC,2)=CSXR(IREC1,1)
CSYR(IREC,2)=CSYR(IREC1,1)
CSZR(IREC,2)=CSZR(IREC1,1)
SNXR(IREC,2)=SNXR(IREC1,1)
TAUPSR(IREC,2)=TAUPSR(IREC1,1)
CPSIR(IREC,2)=CPSIR(IREC1,1)
JJR(IREC,2)=JJR(IREC1,1)
KO(IREC,JJR(IREC,2),2)=KO(IREC1,JJR(IREC1,1),1)
INOUT(IREC,2)=INOUT(IREC1,1)
ENDDO
C
NREC2=NREC2+NREC1
MAXA=MAX0(MAXA,NREC2)
NALL=NALL+NREC2
NREC1=0
IF(NREC2.GT.2000) GO TO 8885
GO TO 8886
8885 WRITE(6,8887)
8887 FORMAT(1X,'ERROR : DIMENSION IN THE RECOIL LOOP ',
& 'MUST BE INCREASED')
8886 CONTINUE
C
C PROCESS THE PARTICLES IN LIST 2
C
C FIND LAYER
C
DO IREC1=1,NREC2
LRR(IREC1,2)=MIN0(ISRCHFGT(L,XX(1),1,XR(IREC1,2)),L)
ENDDO
C
C MOVE PARTICLES
C
DO IREC1=1,NREC2
XR(IREC1,2)=XR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2))
& *CSXR(IREC1,2)
YR(IREC1,2)=YR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2))
& *CSYR(IREC1,2)
ZR(IREC1,2)=ZR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2))
& *CSZR(IREC1,2)
CXR(IREC1)=CSXR(IREC1,2)
CYR(IREC1)=CSYR(IREC1,2)
CZR(IREC1)=CSZR(IREC1,2)
SXR(IREC1)=SNXR(IREC1,2)
DEERS(IREC1)=0.D0
TS(IREC1)=0.D0
ENDDO
KK2=KK0R
DO KKR=KK2,0,-1
C
C CHOICE OF COLLISION PARTNERS
C
DO IREC1=1,NREC2
call ranlux(random, 1)
JJR(IREC1,1) = ISRCHFGE(NJ(LRR(IREC1,2)),COM(1,LRR(IREC1,2))
& ,1,DBLE(random(1)))+JT(LRR(IREC1,2))
ENDDO
DO IREC1=1,NREC2
call ranlux(ran2, 2)
PHIPR=PI2*DBLE(ran2(1))
CPHIR(IREC1,2)=DCOS(PHIPR)
SPHIR(IREC1,2)=DSIN(PHIPR)
PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(ran2(2))+KKR)
X2(IREC1)=XR(IREC1,2)-PR(IREC1)*CPHIR(IREC1,2)*SXR(IREC1)
PR(IREC1)=CVMGT(1.D10,PR(IREC1),X2(IREC1).LT.0.0D0 .OR
& .X2(IREC1).GT.TT)
BR(IREC1)=PR(IREC1)/A(JJR(IREC1,2),JJR(IREC1,1))
EPSR(IREC1)=F(JJR(IREC1,2),JJR(IREC1,1))*ER(IREC1,2)
ENDDO
CALL SCOPY(NREC2,BR,1,RR,1)
IVMIN=1
IVMAX=NREC2
C
C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.)
C
IF(IPOTR.NE.POT_KRC) GO TO 4201
C KR-C POTENTIAL
C CALL MAGICKRC(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2)
205 DO IV=IVMIN,IVMAX
IF(RR(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 205 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
RR(IV)=0.00001D0
ENDIF
EX1R(IV)=DEXP(-.278544D0*RR(IV))
EX2R(IV)=DEXP(-.637174D0*RR(IV))
EX3R(IV)=DEXP(-1.919249D0*RR(IV))
RRR1=1./RR(IV)
VR(IV)=(.190945D0*EX1R(IV)+.473674D0*EX2R(IV)+.335381D0
& *EX3R(IV))*RRR1
FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV)
V1R(IV)=-(VR(IV)+.053186584080D0*EX1R(IV)+.301812757276D0
& *EX2R(IV)+.643679648869D0*EX3R(IV))*RRR1
FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/EPSR(IV
& )-1.D0
Q=FR/FR1
RR(IV)=RR(IV)-Q
TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0
ENDDO
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 207
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 207
GO TO 205
207 DO IV=1,NREC2
ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV))
SQE=DSQRT(EPSR(IV))
CC=(.235800D0+SQE)/(.126000D0+SQE)
AA=2.D0*EPSR(IV)*(1.D0+(1.0144D0/SQE))*BR(IV)**CC
FF=(DSQRT(AA*AA+1.D0)-AA)*((69350.D0+EPSR(IV)) /(83550.D0
& +EPSR(IV)))
DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0)
C2R(IV)=DMIN1(1.0D0,C*C)
S2R(IV)=1.D0-C2R(IV)
ENDDO
GO TO 4203
4201 IF(IPOTR.NE.POT_MOLIERE) GO TO 4202
C MOLIERE POTENTIAL
C CALL MAGICMOL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2)
4205 DO IV=IVMIN,IVMAX
IF(RR(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 4205 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
RR(IV)=0.00001D0
ENDIF
EX1R(IV)=DEXP(-.3D0*RR(IV))
EX2R(IV)=DEXP(-1.2D0*RR(IV))
EX3R(IV)=DEXP(-6.0D0*RR(IV))
RRR1=1.D0/RR(IV)
VR(IV)=(.35D0*EX1R(IV)+.55D0*EX2R(IV)+.10D0*EX3R(IV))*RRR1
FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV)
V1R(IV)=-(VR(IV)+.105D0*EX1R(IV)+ .66D0*EX2R(IV)+.6D0
& *EX3R(IV))*RRR1
FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/
& EPSR(IV)-1.D0
Q=FR/FR1
RR(IV)=RR(IV)-Q
TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0
ENDDO
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 4207
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 4207
GO TO 4205
4207 DO IV=1,NREC2
ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV))
SQE=DSQRT(EPSR(IV))
CC=(.009611D0+SQE)/(.005175D0+SQE)
AA=2.D0*EPSR(IV)*(1.D0+(0.6743D0/SQE))*BR(IV)**CC
FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPSR(IV))/(10.+EPSR(IV))
& )
DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0)
C2R(IV)=DMIN1(1.0D0,C*C)
S2R(IV)=1.D0-C2R(IV)
ENDDO
GO TO 4203
4202 IF(IPOTR.NE.POT_ZBL) GO TO 4203
C ZBL POTENTIAL
C CALL MAGICZBL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2)
5205 DO 5206 IV=IVMIN,IVMAX
IF(RR(IV).LT.1.D-20)THEN
WRITE(99,'(A70)'
& )'in DO 5205 R(IV)<1.D-20 -> 0.00001D0 gesetzt'
RR(IV)=0.00001D0
ENDIF
EX1R(IV)=DEXP(-.20162D0*RR(IV))
EX2R(IV)=DEXP(-.40290D0*RR(IV))
EX3R(IV)=DEXP(-.94229D0*RR(IV))
EX4R(IV)=DEXP(-3.1998D0*RR(IV))
RRR1=1./RR(IV)
VR(IV)=(.02817D0*EX1R(IV)+.28022D0*EX2R(IV)+.50986D0*EX3R(IV
& )+.18175D0*EX4R(IV))*RRR1
FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV)
V1R(IV)=-(VR(IV)+.0056796354D0*EX1R(IV)+.112900638D0*EX2R(IV
& )+.4804359794D0*EX3R(IV)+.581563650D0*EX4R(IV))*RRR1
FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/
& EPSR(IV)-1.D0
Q=FR/FR1
RR(IV)=RR(IV)-Q
TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0
5206 CONTINUE
C GET MAX AND MIN INDEX OF TEST FAILURES
IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1)
IF(IVMIN.GT.IVMAX) GO TO 5207
IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1)
IF(IVMIN.GT.IVMAX) GO TO 5207
GO TO 5205
5207 DO IV=1,NREC2
ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV))
SQE=DSQRT(EPSR(IV))
CC=(.011615D0+SQE)/(.0071222D0+SQE)
AA=2.*EPSR(IV)*(1.D0+(0.99229D0/SQE))*BR(IV)**CC
FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPSR(IV)) /(14.813D0
& +EPSR(IV)))
DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0)
C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0)
C2R(IV)=DMIN1(1.0D0,C*C)
S2R(IV)=1.D0-C2R(IV)
ENDDO
4203 CONTINUE
C
DO IREC1=1,NREC2
T(IREC1)=ER(IREC1,2)*S2R(IREC1)*EC(JJR(IREC1,2),JJR(IREC1,1)
& )
TS(IREC1)=TS(IREC1)+T(IREC1)
T1=CVMGT(T(IREC1),0.D0,KKR.EQ.3)
TR1=TR1+T1
DEEORR=CVMGT(0.D0,KOR(JJR(IREC1,2),JJR(IREC1,1))
& * DSQRT(ER(IREC1,2))*EX1R(IREC1),KDEE2.EQ.STOP_LS)
DEERS(IREC1)=DEERS(IREC1)+DEEORR
TAUR(IREC1)=CVMGT(PR(IREC1)*DSQRT(S2R(IREC1)/C2R(IREC1)),0
& .D0,KKR.EQ.0)
TAUR(IREC1)=DMIN1(TAUR(IREC1),LM(LRR(IREC1,2)))
CTR(IREC1)=C2R(IREC1)+C2R(IREC1)-1.D0
STR(IREC1)=DSQRT(1.D0-CTR(IREC1)*CTR(IREC1))
CUR(IREC1) = CTR(IREC1)+MU(JJR(IREC1,2),JJR(IREC1,1))
CUR(IREC1) = CVMGT(CUR(IREC1),1.0D-8,DABS(CUR(IREC1)).GE.1
& .0D-8)
TAR=STR(IREC1)/CUR(IREC1)
TAR2=1./DSQRT(1.D0+TAR*TAR)
CPSIR(IREC1,2)=TAR2
SPSIR(IREC1,2)=TAR*TAR2
ENDDO
CALL DIRCOS(CSXR(1,2),CSYR(1,2),CSZR(1,2),SNXR(1,2), CPSIR(1,2)
& ,SPSIR(1,2),CPHIR(1,2),SPHIR(1,2),NREC2)
ENDDO
C
C CREATE SECONDARY KNOCK-ON ATOMS
C
DO 246 IREC1=1,NREC2
IF(T(IREC1).LE.ERC) GO TO 246
IF(X2(IREC1).GT.RD.AND.X2(IREC1).LT.RT) GO TO 246
NREC1=NREC1+1
ER(NREC1,1)=T(IREC1)-EP(JJR(IREC1,1))
XR(NREC1,1)=X2(IREC1)
YR(NREC1,1)=YR(IREC1,2)-PR(IREC1)*(SPHIR(IREC1,2)*CZR(IREC1)-
& CPHIR(IREC1,2)*CYR(IREC1)*CXR(IREC1))/SXR(IREC1)
ZR(NREC1,1)=ZR(IREC1,2)+PR(IREC1)*(SPHIR(IREC1,2)*CYR(IREC1)+
& CPHIR(IREC1,2)*CXR(IREC1)*CZR(IREC1))/SXR(IREC1)
CSXR(NREC1,1)=CXR(IREC1)
CSYR(NREC1,1)=CYR(IREC1)
CSZR(NREC1,1)=CZR(IREC1)
SNXR(NREC1,1)=SXR(IREC1)
CPHIR(NREC1,1)=-CPHIR(IREC1,2)
SPHIR(NREC1,1)=-SPHIR(IREC1,2)
CTR(NREC1)=DMIN1(CTR(IREC1),.99999999D0)
TAR=STR(IREC1)/(1.D0-CTR(NREC1))
TAR2=1./DSQRT(1.D0+TAR*TAR)
CPSIR(NREC1,1)=TAR2
SPSIR(NREC1,1)=TAR*TAR2
TAUPSR(NREC1,1)=0.D0
KO(NREC1,JJR(IREC1,1),1)=0
INOUT(NREC1,1)=INOUT(IREC1,2)
JJR(NREC1,1)=JJR(IREC1,1)
NSA=NSA+1
246 CONTINUE
C
C INELASTIC ENERGY LOSS
C
DO IREC1=1,NREC2
ASIGTR(IREC1)=(LM(LRR(IREC1,2))-TAUR(IREC1)+ TAUPSR(IREC1,2))
& *ARHO(LRR(IREC1,2))
TAUPSR(IREC1,2)=TAUR(IREC1)*DABS(CPSIR(IREC1,2))
ENDDO
C KDEE2 uses STOP_LS, STOP_OR, and STOP_MIXED in this order.
GO TO(115,116,117),KDEE2
115 DO IREC1=1,NREC2
DEER(IREC1)=CVMGT(0.D0,KLM(LRR(IREC1,2), JJR(IREC1,2))
& *ASIGTR(IREC1)*DSQRT(ER(IREC1,2)), XR(IREC1,2).LT.HLM.OR
& .XR(IREC1,2).GT.HLMT)
ENDDO
GO TO 242
116 DO IREC1=1,NREC2
DEER(IREC1)=DEERS(IREC1)
ENDDO
GO TO 242
117 DO IREC1=1,NREC2
DEER(IREC1)=CVMGT(DEERS(IREC1),.5*(KLM(LRR(IREC1,2),JJR(IREC1,2
& ))*ASIGTR(IREC1)*DSQRT(ER(IREC1,2))+DEERS(IREC1)),XR(IREC1
& ,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT)
ENDDO
242 CONTINUE
C
DO IREC1=1,NREC2
DELR=DMAX1(1.0D-20,TS(IREC1)+DEER(IREC1))
TS(IREC1)=CVMGT(ER(IREC1,2)*TS(IREC1)/DELR,TS(IREC1) ,DELR.GT
& .ER(IREC1,2))
DEER(IREC1)=CVMGT(ER(IREC1,2)*DEER(IREC1)/DELR,DEER(IREC1)
& ,DELR.GT.ER(IREC1,2))
ENDDO
DO 252 IREC1=1,NREC2
I=MAX0(MIN0(IDINT(X2(IREC1)/CW+1.D0),100),1)
DENTR(I)=DENTR(I)+TS(IREC1)
DMGNR(I)=DMGNR(I)+T(IREC1)
IONR(I)=IONR(I)+DEER(IREC1)
IF(T(IREC1).LE.DI(JJR(IREC1,1))) GO TO 84
ELGDR(I)=ELGDR(I)+T(IREC1)
ICDR(I,JJR(IREC1,2))=ICDR(I,JJR(IREC1,2))+1
ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))= ICDIRI(I,JJR(IREC1,2)
& ,JJR(IREC1,1))+1
GO TO 252
84 PHONR(I)=PHONR(I)+T(IREC1)
252 CONTINUE
DO IREC1=1,NREC2
ICDIR=ICDIR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1)))
& )
ICSBR=ICSBR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.SB(1)))
ICSUMR=ICSUMR+IDINT(CVMGT(1.D0,0.D0,TS(IREC1).GT.0.D0))
ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= ICDIRJ(JJR(IREC1,2)
& ,JJR(IREC1,1)) +IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT
& .DI(JJR(IREC1,1))))
TEL=TEL+TS(IREC1)
TR=TR+T(IREC1)
EI=EI+DEER(IREC1)
ER(IREC1,2)=ER(IREC1,2)-DEER(IREC1)-TS(IREC1)+1.0D-10
XR(IREC1,2)=XR(IREC1,2)-TAUR(IREC1)*CXR(IREC1)
YR(IREC1,2)=YR(IREC1,2)-TAUR(IREC1)*CYR(IREC1)
ZR(IREC1,2)=ZR(IREC1,2)-TAUR(IREC1)*CZR(IREC1)
TESTR(IREC1)=ER(IREC1,2).LE.SFE.OR.XR(IREC1,2).LT.(-SU) .OR
& .XR(IREC1,2).GT.SUT .OR.(XR(IREC1,2).GT.RD.AND.XR(IREC1,2)
& .LT.RT)
ENDDO
C
C CHECK TO SEE IF ANY RECOIL ATOM IS SPUTTERED OR
C IF THE ENERGY OF ANY RECOIL ATOM IS TOO LOW
C
IVMIN=1+ILLZ(NREC2,TESTR,1)
IF(IVMIN.GT.NREC2) GO TO 247
IVMAX=NREC2-ILLZ(NREC2,TESTR,-1)
C
DO 248 IREC1=IVMIN,IVMAX
249 IF(IREC1.GT.NREC2) GO TO 247
IF(XR(IREC1,2).LT.(-SU)) GO TO 251
IF(XR(IREC1,2).GT.SUT) GO TO 250
IF(ER(IREC1,2).LE.ERC) GO TO 255
IF(XR(IREC1,2).GT.RD.AND.XR(IREC1,2).LT.RT) GO TO 255
GO TO 248
251 ENOR=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1)
IF(ENOR.GT.SB(1)) GO TO 254
C
C RECOIL ATOM IS REFLECTED BACK INTO THE SOLID BY THE
C POTENTIAL BARRIER
C
XR(IREC1,2)=-1.D0*SU
CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2)
KIS=KIS+1
GO TO 248
C
C RECOIL ATOM IS SPUTTERED (BACKWARD)
C
254 ESP=ER(IREC1,2)-SB(1)
CC254 ESP=ER(IREC1,2)-SB(LRR(IREC1,2))
C
C NUMBER, ENERGY AND MOMENTS OF ALL SPUTTERED PARTICLES
C
IBSP(JJR(IREC1,2))=IBSP(JJR(IREC1,2))+1
EBSP(JJR(IREC1,2))=EBSP(JJR(IREC1,2))+ESP
SPE2=ESP*ESP
SPE3=SPE2*ESP
SPE2S(JJR(IREC1,2))=SPE2S(JJR(IREC1,2))+SPE2
SPE3S(JJR(IREC1,2))=SPE3S(JJR(IREC1,2))+SPE3
SPE4S(JJR(IREC1,2))=SPE4S(JJR(IREC1,2))+SPE2*SPE2
SPE5S(JJR(IREC1,2))=SPE5S(JJR(IREC1,2))+SPE3*SPE2
SPE6S(JJR(IREC1,2))=SPE6S(JJR(IREC1,2))+SPE3*SPE3
IF(ESP.LT.0.1) GO TO 256
IBSPL(JJR(IREC1,2))=IBSPL(JJR(IREC1,2))+1
SPE1SL(JJR(IREC1,2))=SPE1SL(JJR(IREC1,2))+DLOG10(DABS(ESP))
SPE2L=(DLOG10(DABS(ESP)))**2.D0
SPE3L=SPE2L*DLOG10(DABS(ESP))
SPE2SL(JJR(IREC1,2))=SPE2SL(JJR(IREC1,2))+SPE2L
SPE3SL(JJR(IREC1,2))=SPE3SL(JJR(IREC1,2))+SPE3L
SPE4SL(JJR(IREC1,2))=SPE4SL(JJR(IREC1,2))+SPE2L*SPE2L
SPE5SL(JJR(IREC1,2))=SPE5SL(JJR(IREC1,2))+SPE3L*SPE2L
SPE6SL(JJR(IREC1,2))=SPE6SL(JJR(IREC1,2))+SPE3L*SPE3L
256 CONTINUE
C
C SURFACE REFRACTION
C
EXIR=DSQRT((ENOR-SB(1))/ESP)
IF ( EXIR .GE. 1.D0 ) EXIR = .999999D0
C
C TOTAL ANGULAR DISTRIBUTIONS
C
IAG=IDINT(EXIR*20.D0+1.D0)
KADS(IAG)=KADS(IAG)+1
KADSJ(IAG,JJR(IREC1,2))=KADSJ(IAG,JJR(IREC1,2))+1
C
C 4 GROUPS:ION IN , PKA ;ION IN , SKA ;ION OUT, PKA ;ION OUT, SKA
C
KOI=KO(IREC1,JJR(IREC1,2),2)
IF(INOUT(IREC1,2).EQ.-1) GO TO 61
IF(KOI.EQ.0) GO TO 62
ISPIP(JJR(IREC1,2))=ISPIP(JJR(IREC1,2))+1
ESPIP(JJR(IREC1,2))=ESPIP(JJR(IREC1,2))+ESP
KADRIP(IAG,JJR(IREC1,2))=KADRIP(IAG,JJR(IREC1,2))+1
GO TO 164
62 ISPIS(JJR(IREC1,2))=ISPIS(JJR(IREC1,2))+1
ESPIS(JJR(IREC1,2))=ESPIS(JJR(IREC1,2))+ESP
KADRIS(IAG,JJR(IREC1,2))=KADRIS(IAG,JJR(IREC1,2))+1
GO TO 164
61 IF(KOI.EQ.0) GO TO 163
ISPOP(JJR(IREC1,2))=ISPOP(JJR(IREC1,2))+1
ESPOP(JJR(IREC1,2))=ESPOP(JJR(IREC1,2))+ESP
KADROP(IAG,JJR(IREC1,2))=KADROP(IAG,JJR(IREC1,2))+1
GO TO 164
163 ISPOS(JJR(IREC1,2))=ISPOS(JJR(IREC1,2))+1
ESPOS(JJR(IREC1,2))=ESPOS(JJR(IREC1,2))+ESP
KADROS(IAG,JJR(IREC1,2))=KADROS(IAG,JJR(IREC1,2))+1
164 CONTINUE
C
C OUTPUT MATRICES OF BACKWARD SPUTTERED ATOMS
C
IF(JJR(IREC1,2).GT.JT(3)) GO TO 255
IESP = IDINT(ESP*E0DE+2.0D0)
IESP = MIN0(101,IESP)
IESLOG=2
IF(ESP.LE.0.1D0) GO TO 75
IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESP))+3.D0)
IESLOG=MIN0(IESLOG,75)
75 CONTINUE
IF(EXIR.GT.1.0D0)WRITE(99,'(A50)')' EXIR nach Label 75'
IA=IDINT(DAW*DACOS(EXIR)+2.D0)
IAGS = IAG+1
IG=2
SXR(IREC1)=DMAX1(SXR(IREC1),1.0D-12)
U=CYR(IREC1)/SXR(IREC1)
IF(DABS(U).GT.1.D0) U = SIGN(1.D0,U)
IF(U.GT.1.0D0)WRITE(99,'(A50)')' U vor Label 182 CON..'
ACS=DACOS(U)
IG=IDINT(DGW*ACS+2.D0)
IGG=IDINT(DGIK*ACS+1.D0)
IGG=MAX0(MIN0(NGIK,IGG),1)
MAGS(IG,IAGS,JJR(IREC1,2))=MAGS(IG,IAGS,JJR(IREC1,2))+1
MAGS(IG,22,JJR(IREC1,2))=MAGS(IG,22,JJR(IREC1,2))+1
MAGS(NG,IAGS,JJR(IREC1,2))=MAGS(NG,IAGS,JJR(IREC1,2))+1
MAGS(NG,22,JJR(IREC1,2))=MAGS(NG,22,JJR(IREC1,2))+1
MEAS(IESP,IAGS,JJR(IREC1,2))=MEAS(IESP,IAGS,JJR(IREC1,2))+1
MEAS(102,IAGS,JJR(IREC1,2))=MEAS(102,IAGS,JJR(IREC1,2))+1
MEAS(IESP,22,JJR(IREC1,2))=MEAS(IESP,22,JJR(IREC1,2))+1
MEAS(102,22,JJR(IREC1,2))=MEAS(102,22,JJR(IREC1,2))+1
MEASL(IESLOG,IAG,JJR(IREC1,2))=MEASL(IESLOG,IAG,JJR(IREC1,2))+1
MEASL(IESLOG,21,JJR(IREC1,2))=MEASL(IESLOG,21,JJR(IREC1,2))+1
MEASL(75,IAG,JJR(IREC1,2))=MEASL(75,IAG,JJR(IREC1,2))+1
MEASL(75,21,JJR(IREC1,2))=MEASL(75,21,JJR(IREC1,2))+1
IF(ALPHA.LT.1.0) GO TO 181
MEAGS(IESP,IGG,IAGS,JJR(IREC1,2)) = MEAGS(IESP,IGG,IAGS
& ,JJR(IREC1,2))+1
MEAGS(102,IGG,IAGS,JJR(IREC1,2)) = MEAGS(102,IGG,IAGS
& ,JJR(IREC1,2))+1
MEAGS(IESP,IGG,22,JJR(IREC1,2)) = MEAGS(IESP,IGG,22,JJR(IREC1
& ,2))+1
MEAGS(102,IGG,22,JJR(IREC1,2)) = MEAGS(102,IGG,22,JJR(IREC1
& ,2))+1
MAGSA(IG,IA,JJR(IREC1,2)) = MAGSA(IG,IA,JJR(IREC1,2))+1
MAGSA(IG,32,JJR(IREC1,2)) = MAGSA(IG,32,JJR(IREC1,2))+1
MAGSA(NG,IA,JJR(IREC1,2)) = MAGSA(NG,IA,JJR(IREC1,2))+1
MAGSA(NG,32,JJR(IREC1,2)) = MAGSA(NG,32,JJR(IREC1,2))+1
181 CONTINUE
GO TO 255
C
250 ENORT=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1)
IF(ENORT.GT.SB(L)) GO TO 257
C
C RECOIL ATOM IS REFLECTED BACK INTO THE SOLID BY THE
C POTENTIAL BARRIER
C
XR(IREC1,2)=SUT
CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2)
KIST=KIST+1
GO TO 248
C
C RECOIL ATOM IS SPUTTERED (TRANSMISSION)
C
257 ESPT=ER(IREC1,2)-SB(L)
C
C NUMBER AND ENERGY OF ALL SPUTTERED PARTICLES
C
ITSP(JJR(IREC1,2))=ITSP(JJR(IREC1,2))+1
ETSP(JJR(IREC1,2))=ETSP(JJR(IREC1,2))+ESPT
C
C SURFACE REFRACTION
C
EXIRT=DSQRT((ENORT-SB(L))/ESPT)
IF ( EXIRT .GE. 1.D0 ) EXIRT = .999999D0
C
C TOTAL ANGULAR DISTRIBUTIONS
C
IAG=IDINT(EXIRT*20.D0+1.D0)
KADST(IAG)=KADST(IAG)+1
KDSTJ(IAG,JJR(IREC1,2))=KDSTJ(IAG,JJR(IREC1,2))+1
C BUGFIX: transmission-sputter matrices below use IG.
C Compute it here instead of reusing a stale value.
IG=2
SXR(IREC1)=DMAX1(SXR(IREC1),1.0D-12)
U=CYR(IREC1)/SXR(IREC1)
IF(DABS(U).GT.1.D0) U = SIGN(1.D0,U)
ACS=DACOS(U)
IG=IDINT(DGW*ACS+2.D0)
C
C 4 GROUPS:ION IN , PKA ;ION IN , SKA ;ION OUT, PKA ;ION OUT, SKA
C
KOI=KO(IREC1,JJR(IREC1,2),2)
IF(INOUT(IREC1,2).EQ.-1) GO TO 85
IF(KOI.EQ.0) GO TO 86
ISPIPT(JJR(IREC1,2))=ISPIPT(JJR(IREC1,2))+1
ESPIPT(JJR(IREC1,2))=ESPIPT(JJR(IREC1,2))+ESPT
GO TO 88
86 ISPIST(JJR(IREC1,2))=ISPIST(JJR(IREC1,2))+1
ESPIST(JJR(IREC1,2))=ESPIST(JJR(IREC1,2))+ESPT
GO TO 88
85 IF(KOI.EQ.0) GO TO 87
ISPOPT(JJR(IREC1,2))=ISPOPT(JJR(IREC1,2))+1
ESPOPT(JJR(IREC1,2))=ESPOPT(JJR(IREC1,2))+ESPT
GO TO 88
87 ISPOST(JJR(IREC1,2))=ISPOST(JJR(IREC1,2))+1
ESPOST(JJR(IREC1,2))=ESPOST(JJR(IREC1,2))+ESPT
88 CONTINUE
C
C OUTPUT MATRICES OF FORWARD SPUTTERED ATOMS
C
JRT=JJR(IREC1,2)
IF(L.EQ.3) JRT=JJR(IREC1,2)-NJ(1)
IF(JRT.LT.1) GO TO 255
IESPT = IDINT(ESPT*E0DE+2.0D0)
IESPT = MIN0(101,IESPT)
IESLOG=2
IF(ESPT.LE.0.1D0) GO TO 76
IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESPT))+3.D0)
IESLOG=MIN0(IESLOG,75)
76 CONTINUE
IAGS = IAG+1
MAGST(IG,IAGS,JRT) = MAGST(IG,IAGS,JRT)+1
MAGST(IG,22,JRT) = MAGST(IG,22,JRT)+1
MAGST(NG,IAGS,JRT) = MAGST(NG,IAGS,JRT)+1
MAGST(NG,22,JRT) = MAGST(NG,22,JRT)+1
MEAST(IESPT,IAGS,JRT) = MEAST(IESPT,IAGS,JRT)+1
MEAST(102,IAGS,JRT) = MEAST(102,IAGS,JRT)+1
MEAST(IESPT,22,JRT) = MEAST(IESPT,22,JRT)+1
MEAST(102,22,JRT) = MEAST(102,22,JRT)+1
MEASTL(IESLOG,IAG,JRT) = MEASTL(IESLOG,IAG,JRT)+1
MEASTL(IESLOG,21,JRT) = MEASTL(IESLOG,21,JRT)+1
MEASTL(75,IAG,JRT) = MEASTL(75,IAG,JRT)+1
MEASTL(75,21,JRT) = MEASTL(75,21,JRT)+1
C
C REARRANGEMENT OF PARTICLES IN LIST 2
C
255 ER(IREC1,2)=ER(NREC2,2)
XR(IREC1,2)=XR(NREC2,2)
YR(IREC1,2)=YR(NREC2,2)
ZR(IREC1,2)=ZR(NREC2,2)
CSXR(IREC1,2)=CSXR(NREC2,2)
CSYR(IREC1,2)=CSYR(NREC2,2)
CSZR(IREC1,2)=CSZR(NREC2,2)
SNXR(IREC1,2)=SNXR(NREC2,2)
CPHIR(IREC1,2)=CPHIR(NREC2,2)
SPHIR(IREC1,2)=SPHIR(NREC2,2)
CPSIR(IREC1,2)=CPSIR(NREC2,2)
SPSIR(IREC1,2)=SPSIR(NREC2,2)
TAUPSR(IREC1,2)=TAUPSR(NREC2,2)
JJR(IREC1,2)=JJR(NREC2,2)
KO(IREC1,JJR(IREC1,2),2)=KO(NREC2,JJR(NREC2,2),2)
INOUT(IREC1,2)=INOUT(NREC2,2)
NREC2=NREC2-1
C
IF(IREC1.EQ.NREC2+1) GO TO 247
C THE NREC2 PARTICLE FAILS THE TEST
IF(NREC2+1.GT.IVMAX) GO TO 248
GO TO 249
248 CONTINUE
247 CONTINUE
IF(NREC1+NREC2.EQ.0) GO TO 27
IF(NREC2.GE.NUM.OR.IH1.EQ.0) GO TO 83
C
C END OF RECOIL ATOM SECTION
C
27 CONTINUE
IF(IH1.EQ.0.AND.IH.EQ.NH) GO TO 140
C
C=======================================================================
C PROJECTILE STOP, BACKSCATTER AND TRANSMISSION TESTS
C=======================================================================
C PROJECTILE CANDIDATE FOR REFLECTION
C
DO IV=1,IH1
E(IV)=E(IV)-DEE(IV)-DENS(IV)+1.0D-10
X(IV)=X(IV)-TAU(IV)*CX(IV)
Y(IV)=Y(IV)-TAU(IV)*CY(IV)
Z(IV)=Z(IV)-TAU(IV)*CZ(IV)
PL(IV)=PL(IV)-TAU(IV)
TEST(IV)=E(IV).LE.EF.OR.X(IV).LT.-1.D0*SU.OR.X(IV).GT.SUT
ENDDO
IVMIN=1+ILLZ(IH1,TEST,1)
IF(IVMIN.GT.IH1) THEN
C BUGFIX: no projectile needs special handling; ensure IVMAX
C is defined before the shared label-90 bookkeeping.
IVMAX=0
GO TO 90
ENDIF
IVMAX=IH1-ILLZ(IH1,TEST,-1)
DO 120 IV=IVMIN,IVMAX
160 IF(IV.GT.IH1) GO TO 90
IF(X(IV).LT.-SU) GO TO 8
IF(X(IV).GT.SUT) GO TO 9
IF(E(IV).GT.EF) GO TO 125
IF(E(IV).GT.ESB.AND.X(IV).LT.0.D0) GO TO 125
IF(E(IV).GT.ESB.AND.X(IV).GT.TT) GO TO 125
C
C PROJECTILE HAS STOPPED (PATHLENGTH,RANGE,SPREAD AND MOMENTS)
C
IP = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), MAXD), 1)
IPL(IP)=IPL(IP)+1
I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), MAXD1), 0)
IRP(I1)=IRP(I1)+1
C
C Calculate the stopped particles in each layer
C
LowTiefe = 0.D0
UpTiefe = DX(1)
DO laufzahl=1,l
IF(X(IV).GT.LowTiefe.AND.X(IV).LE.UpTiefe) THEN
Number_in_Layer(laufzahl)=Number_in_Layer(laufzahl)+1
ENDIF
LowTiefe = UpTiefe
UpTiefe = UpTiefe+DX(laufzahl+1)
ENDDO
PL2=PL(IV)*PL(IV)
PL3=PL2*PL(IV)
PLSUM=PLSUM+PL(IV)
PL2SUM=PL2SUM+PL2
PL3SUM=PL3SUM+PL3
PL4SUM=PL4SUM+PL2*PL2
PL5SUM=PL5SUM+PL3*PL2
PL6SUM=PL6SUM+PL3*PL3
IF(X(IV).LT.0.D0.OR.X(IV).GT.TT) GO TO 111
XQ=X(IV)*X(IV)
XQ3=XQ*X(IV)
XSUM=XSUM+X(IV)
X2SUM=X2SUM+XQ
X3SUM=X3SUM+XQ3
X4SUM=X4SUM+XQ*XQ
X5SUM=X5SUM+XQ3*XQ
X6SUM=X6SUM+XQ3*XQ3
RQ=Y(IV)*Y(IV)+Z(IV)*Z(IV)
RQW=DSQRT(RQ)
RQ3=RQ*RQW
RSUM=RSUM+RQW
R2SUM=R2SUM+RQ
R3SUM=R3SUM+RQ3
R4SUM=R4SUM+RQ*RQ
R5SUM=R5SUM+RQ3*RQ
R6SUM=R6SUM+RQ3*RQ3
111 CONTINUE
ENUCLI=ENUCLI+ENUCL(IV)
ENL2I=ENL2I+ENUCL(IV)*ENUCL(IV)
ENUCL(IV)=0.D0
EINELI=EINELI+EINEL(IV)
EIL2I=EIL2I+EINEL(IV)*EINEL(IV)
EINEL(IV)=0.D0
GO TO 110
8 ENO=E(IV)*CX(IV)*CX(IV)
IF(ENO.LE.ESB) GO TO 24
C
C PROJECTILE IS BACKSCATTERED
C
IB=IB+1
ES=E(IV)-ESB
ESQ=ES*ES
ES3=ESQ*ES
EB=EB+ES
EB2SUM=EB2SUM+ESQ
EB3SUM=EB3SUM+ES3
EB4SUM=EB4SUM+ESQ*ESQ
EB5SUM=EB5SUM+ES3*ESQ
EB6SUM=EB6SUM+ES3*ES3
IF(ES.LT.0.1D0) GO TO 112
IBL=IBL+1
ESQL=(DLOG10(DABS(ES)))**2.D0
ES3L=ESQL*DLOG10(DABS(ES))
EB1SUL=EB1SUL+DLOG10(DABS(ES))
EB2SUL=EB2SUL+ESQL
EB3SUL=EB3SUL+ES3L
EB4SUL=EB4SUL+ESQL*ESQL
EB5SUL=EB5SUL+ES3L*ESQL
EB6SUL=EB6SUL+ES3L*ES3L
112 CONTINUE
IPB = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1)
IPLB(IPB)=IPLB(IPB)+1
PLQB=PL(IV)*PL(IV)
PL3B=PLQB*PL(IV)
PLSB=PLSB+PL(IV)
PL2SB=PL2SB+PLQB
PL3SB=PL3SB+PL3B
PL4SB=PL4SB+PLQB*PLQB
PL5SB=PL5SB+PL3B*PLQB
PL6SB=PL6SB+PL3B*PL3B
ENUCLB=ENUCLB+ENUCL(IV)
ENL2B=ENL2B+ENUCL(IV)*ENUCL(IV)
ENUCL(IV)=0.D0
EINELB=EINELB+EINEL(IV)
EIL2B=EIL2B+EINEL(IV)*EINEL(IV)
EINEL(IV)=0.D0
C
C SURFACE REFRACTION
C
EXI=DSQRT((ENO-ESB)/ES)
exi1s=exi1s+exi
exiq=exi*exi
exic=exiq*exi
exi2s=exi2s+exiq
exi3s=exi3s+exic
exi4s=exi4s+exiq*exiq
exi5s=exi5s+exic*exiq
exi6s=exi6s+exic*exic
C
C DIVISIONS FOR VECTORS AND MATRICES
C
IE = IDINT(E0DE*ES+2.D0)
IE = MAX0( MIN0( IE,NE1), 2)
IERLOG = 2
IF(ES.LE.0.1D0) GO TO 4
IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*ES))+3.D0)
IERLOG=MIN0(IERLOG,75)
4 CONTINUE
IAG=IDINT(EXI*20.D0+1.D0)
IAG = MIN0( IAG, 20)
IAGB = IAG+1
KADB(IAG)=KADB(IAG)+1
IG=2
COSSIN=CY(IV)/SX(IV)
COSSIN=DMIN1(COSSIN,1.D0)
COSSIN=DMAX1(COSSIN,-1.D0)
coss1s=coss1s+cossin
cossq=cossin*cossin
cosst=cossq*cossin
coss2s=coss2s+cossq
coss3s=coss3s+cosst
coss4s=coss4s+cossq*cossq
coss5s=coss5s+cosst*cossq
coss6s=coss6s+cosst*cosst
IF(COSSIN.GT.1.0D0)WRITE(99,'(A50)')' nach coss6s'
AC=DACOS(COSSIN)
IG=IDINT(DAW*AC+2.D0)
IGG=IDINT(DGIK*AC+1.D0)
IF(IGG.GT.NGIK) IGG=NGIK
IPB1=IPB+1
MEABL(IERLOG,IAG) = MEABL(IERLOG,IAG)+1
MEABL(IERLOG,21) = MEABL(IERLOG,21)+1
MEABL(75,IAG) = MEABL(75,IAG)+1
MAGB(IG,IAGB) = MAGB(IG,IAGB)+1
MAGB(NG,IAGB) = MAGB(NG,IAGB)+1
MAGB(IG,22) = MAGB(IG,22)+1
MEAB(IE,IAGB) = MEAB(IE,IAGB)+1
MEAB(NE,IAGB) = MEAB(NE,IAGB)+1
MEAB(IE,22) = MEAB(IE,22)+1
MEAGB(IE,IGG,IAGB) = MEAGB(IE,IGG,IAGB)+1
MEAGB(102,IGG,IAGB) = MEAGB(102,IGG,IAGB)+1
MEAGB(IE,IGG,22) = MEAGB(IE,IGG,22)+1
MEAGB(102,IGG,22) = MEAGB(102,IGG,22)+1
MEPB(IE,IPB1) = MEPB(IE,IPB1)+1
MEPB(NE,IPB1) = MEPB(NE,IPB1)+1
MEPB(IE,102) = MEPB(IE,102)+1
EMA(IG,IAGB) = EMA(IG,IAGB)+ES
EMA(IG,22) = EMA(IG,22)+ES
EMA(NG,IAGB) = EMA(NG,IAGB)+ES
GO TO 110
C
C PROJECTILE IS REFLECTED BACK INTO THE TARGET BY THE SURF. BARRIER
C
24 X(IV)=-1.D0*SU
COSX(IV)=-1.D0*COSX(IV)
KIB=KIB+1
GO TO 125
C
C PROJECTILE IS TRANSMITTED
C
9 ENOT=E(IV)*CX(IV)*CX(IV)
IF(ENOT.LE.ESB) GO TO 517
IT=IT+1
EST=E(IV)-ESB
ET=ET+EST
ETQ=EST*EST
ET3=ETQ*EST
ET2SUM=ET2SUM+ETQ
ET3SUM=ET3SUM+ET3
ET4SUM=ET4SUM+ETQ*ETQ
ET5SUM=ET5SUM+ET3*ETQ
ET6SUM=ET6SUM+ET3*ET3
IPT = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1)
C BUGFIX: IPT is computed for transmitted projectiles here.
C IP belongs to different logic and may be uninitialised/stale.
IPLT(IPT)=IPLT(IPT)+1
PLQT=PL(IV)*PL(IV)
PL3T=PLQT*PL(IV)
PLST=PLST+PL(IV)
PL2ST=PL2ST+PLQT
PL3ST=PL3ST+PL3T
PL4ST=PL4ST+PLQT*PLQT
PL5ST=PL5ST+PL3T*PLQT
PL6ST=PL6ST+PL3T*PL3T
ENUCLT=ENUCLT+ENUCL(IV)
ENL2T=ENL2T+ENUCL(IV)*ENUCL(IV)
ENUCL(IV)=0.D0
EINELT=EINELT+EINEL(IV)
EIL2T=EIL2T+EINEL(IV)*EINEL(IV)
EINEL(IV)=0.D0
C
C SURFACE REFRACTION
C
EXI=DSQRT((ENOT-ESB)/EST)
C
C DIVISIONS FOR VECTORS AND MATRICES
C
IE=IDINT(E0DE*EST+2.D0)
IERLOG = 2
IF(EST.LE.0.1D0) GO TO 5
IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*EST))+3.D0)
IERLOG=MIN0(IERLOG,75)
5 CONTINUE
IAG=IDINT(EXI*20.D0+1.D0)
IAG = MIN0( IAG, 20)
IAGB = IAG+1
KADT(IAG)=KADT(IAG)+1
IG=2
COSSIN=CY(IV)/SX(IV)
COSSIN=DMIN1(COSSIN,1.D0)
COSSIN=DMAX1(COSSIN,-1.D0)
IF(COSSIN.GT.1.0D0) WRITE(99,'(A50)')' nach COSSIN'
AC=DACOS(COSSIN)
IG=IDINT(DAW*AC+2.D0)
IGG=IDINT(DGIK*AC+1.D0)
IF(IGG.GT.NGIK) IGG=NGIK
MEATL(IERLOG,IAG) = MEATL(IERLOG,IAG)+1
MEATL(IERLOG,21) = MEATL(IERLOG,21)+1
MEATL(75,IAG) = MEATL(75,IAG)+1
MAGT(IG,IAGB) = MAGT(IG,IAGB)+1
MAGT(NG,IAGB) = MAGT(NG,IAGB)+1
MAGT(IG,22) = MAGT(IG,22)+1
MEAT(IE,IAGB) = MEAT(IE,IAGB)+1
MEAT(NE,IAGB) = MEAT(NE,IAGB)+1
MEAT(IE,22) = MEAT(IE,22)+1
MEAGT(IE,IGG,IAGB) = MEAGT(IE,IGG,IAGB)+1
MEAGT(102,IGG,IAGB) = MEAGT(102,IGG,IAGB)+1
MEAGT(IE,IGG,22) = MEAGT(IE,IGG,22)+1
MEAGT(102,IGG,22) = MEAGT(102,IGG,22)+1
MEPT(IE,IPT) = MEPT(IE,IPT)+1
MEPT(NE,IPT) = MEPT(NE,IPT)+1
MEPT(IE,102) = MEPT(IE,102)+1
C BUGFIX: this is the transmission branch; use EST, not ES
C from the backscattering branch.
EMAT(IG,IAGB) = EMAT(IG,IAGB)+EST
EMAT(IG,22) = EMAT(IG,22)+EST
EMAT(NG,IAGB) = EMAT(NG,IAGB)+EST
GO TO 110
C
C PROJECTILE IS REFLECTED BACK INTO THE TARGET BY THE SURF. BARRIER
C
517 X(IV)=SUT
COSX(IV)=-1.D0*COSX(IV)
KIT=KIT+1
GO TO 125
110 IF(IH.EQ.NH) GO TO 130
IH=IH+1
IF(E0.GE.0.D0) GO TO 702
IF(ALPHA.LT.0.D0) GO TO 703
C
C MAXWELLIAN VELOCITY DISTRIBUTION
C
CALL VELOC(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV))
EMX = EMX+E(IV)
ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0))
me(ne) = me(ne)+1
GO TO 707
C
C MAXWELLIAN ENERGY DISTRIBUTION
C
703 CALL ENERG(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV))
EMX = EMX+E(IV)
GO TO 707
702 IF (EQUAL(Esig,0.D0)) THEN
C FIXED PROJECTILE ENERGY
E(IV)=E0
C GAUSSIAN ENERGY DISTRIBUTION
ELSE
7020 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0)
tryE = tryE+1
IF (Epar.LE.0.D0) THEN
negE = negE+1
GO TO 7020
ENDIF
E(IV)=Epar
ENDIF
TAUPSI(IV)=0.D0
IF(EQUAL(ALPHA,-2.D0)) GO TO 705
IF(EQUAL(ALPHA,-1.D0)) GO TO 706
IF(EQUAL(ALPHASIG,0.D0))THEN
C FIXED PROJECTILE ANGLE
COSX(IV)=CALFA
COSY(IV)=SALFA
COSZ(IV)=0.D0
SINE(IV)=COSY(IV)
ELSE
C
C 1D-GAUSSIAN DISTRIBUTION PROJECTILE ANGLE
C
CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, CALFA
& ,SALFA,BW)
COSX(IV) = CALFA
COSY(IV) = SALFA
COSZ(IV) = 0.D0
SINE(IV) = COSY(IV)
ENDIF
GO TO 707
C
C COSINE ANGLE DISTRIBUTION
C
705 call ranlux(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(1))
COSX(IV)=DSQRT(RTHETA)
SINE(IV)=DSQRT(1.D0-RTHETA)
COSY(IV)=SINE(IV)*DCOS(RPHI)
COSZ(IV)=SINE(IV)*DSIN(RPHI)
GO TO 707
C
C RANDOM DISTRIBUTION
C
706 IF(X0.GT.0.D0) GO TO 709
call ranlux(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(2))
COSX(IV)=1.D0-RTHETA
SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV))
COSY(IV)=SINE(IV)*DSIN(RPHI)
COSZ(IV)=SINE(IV)*DCOS(RPHI)
GO TO 707
709 call ranlux(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(2))
COSX(IV)=1.D0-2.D0*RTHETA
SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV))
COSY(IV)=SINE(IV)*DSIN(RPHI)
COSZ(IV)=SINE(IV)*DCOS(RPHI)
GO TO 708
707 IF(X0.GT.0.D0) GO TO 708
C
C EXTERNAL START
C
SINA=SINE(IV)
COSX(IV)=DSQRT((E(IV)*COSX(IV)*COSX(IV)+ESB)/(E(IV)+ESB))
SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV))
COSY(IV)=COSY(IV)*SINE(IV)/SINA
COSZ(IV)=COSZ(IV)*SINE(IV)/SINA
E(IV)=E(IV)+ESB
C
C LOCUS OF FIRST COLLISION
C
708 LLL(IV)=ISRCHFGT(L,XX(1),1,X0)
call ranlux(random, 1)
RA1=CVMGT(DBLE(random(1)),1.D0,X0.LE.0.D0)
X(IV)=XC+LM(LLL(IV))*RA1*COSX(IV)
Y(IV)=LM(LLL(IV))*RA1*COSY(IV)
Z(IV)=LM(LLL(IV))*RA1*COSZ(IV)
PL(IV)=CVMGT(0.D0,LM(LLL(IV))*RA1,X0.LE.0.D0)
GO TO 120
C
C COUNTING DOWN IH1 , ONLY LESS THAN (NH-IH) HAVE TO BE PROCESSED
C
130 CONTINUE
E(IV)=E(IH1)
COSX(IV)=COSX(IH1)
COSY(IV)=COSY(IH1)
COSZ(IV)=COSZ(IH1)
SINE(IV)=SINE(IH1)
X(IV)=X(IH1)
Y(IV)=Y(IH1)
Z(IV)=Z(IH1)
PL(IV)=PL(IH1)
TAU(IV)=TAU(IH1)
TAUPSI(IV)=TAUPSI(IH1)
CPSI(IV)=CPSI(IH1)
ENUCL(IV)=ENUCL(IH1)
EINEL(IV)=EINEL(IH1)
IH1=IH1-1
IF(IV.EQ.IH1+1) GO TO 90
IF(IH1+1.GT.IVMAX) GO TO 125
GO TO 160
125 CONTINUE
X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV)
Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV)
Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV)
PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0
& .D0)
LLL(IV)=MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L)
120 CONTINUE
90 CONTINUE
C
C INCREMENT OF PROJECTILE ENERGY AND POSITION
C OF PARTICLES NOT HANDLED IN LOOP 120
C
DO 128 IV=1,IVMIN-1
LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L)
128 CONTINUE
DO IV=1,IVMIN-1
X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV)
Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV)
Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV)
PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0
& .D0)
ENDDO
134 CONTINUE
C
C BUGFIX: if no projectile satisfied TEST, IVMIN is IH1+1 and
C IVMAX was not set in this block. All particles were already
C advanced by the loop above, so skip the trailing range.
IF(IVMIN.GT.IH1) GO TO 132
IF(IVMAX.LT.IVMIN) GO TO 132
DO IV=IVMAX+1,IH1
LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L)
X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV)
Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV)
Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV)
PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0
& .D0)
ENDDO
132 CONTINUE
GO TO 1
140 IF(NREC1+NREC2.GT.0) GO TO 83
C
C
C=======================================================================
C WRITE MAIN OUTPUT REPORT
C=======================================================================
C PRINTOUT
C
C
do ima = 5000,1,-1
if(me(ima).ne.0) goto 1010
enddo
ima = 1
1010 ima = MIN0(ima+2,5000)
C I am not sure what is this file for?!
C open(20,file='edist')
C do ne=1,ima
C write(20,1020) ne,me(ne)
C enddo
C 1020 format(1x,2i6)
C close(20)
C
C Calculate the part. reflec. coeff. according to Thomas et al.
C
E0keV=E0/1.D3
EsigkeV=Esig/1.D3
C BUGFIX: single-element target should be decided from NJ(1).
C ZT(1,2) is not read for NJ(1)=1 and may be uninitialised.
IF(NJ(1).EQ.1) THEN
epsilon = 32.55D0*(MT(1,1)/M1)/(1.D0+(MT(1,1)/M1))* 1.D0/(Z1
& *ZT(1,1)*DSQRT(Z1**(2.D0/3.D0)+ZT(1,1)**(2.D0/3.D0)))
& * E0keV
prcoeff = prc(1)*DLOG(prc(2)*epsilon+DEXP(1.D0))/ (1.D0+(PRC(3)
& *epsilon**PRC(4))+(PRC(5)*epsilon**PRC(6)))
ELSE
epsilon = 0.D0
prcoeff = 0.D0
ENDIF
C 2nd CALL DATE_AND_TIME
C
C Write input and beam summary to the main output file.
C
call writeInputSummary(innam,outnam,rgenam,Z1,M1,E0,Esig,
& ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC,NH,RI,RI2,RI3,X0,RD,
& CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT,IPOTR,IRL,L,DX,RHO,
& CK,NJ,ZT,MT,CO,SBE,ED,BE,CH1,CH2,CH3,CH4,CH5,TI,
& ZARG,VELC,EMX,MAXNL,MAXEL,seconds_start_total)
call writeTargetSetupSummary(HLM,HLMT,SU,SUT,XC,RT,SFE,INEL,
& L,LJ,NPROJ,KIB,KIT,MAXA,NALL,NPA,NSA,KIS,KIST,EPS0,Z2,
& M2,ARHO,LM,PDMAX,ASIG,SB,XX,NJ,A1,KOR1,A,KOR,MAXNL,
& MAXNL5)
C
C INTEGRAL IMPLANTATION , SPUTTERING , BACKSCATTERING , TRANSMISSION
C
IIM=NH-IB-IT
YH=DBLE(IIM)
HN=DBLE(NH)
EMV=CVMGT(EMX/HN,E0,E0.LE.0.D0)
EIM=DBLE(IIM)*EMV
DO J=1,LJ
ISPA = ISPA+IBSP(J)
ESPA = ESPA+EBSP(J)
ISPAT = ISPAT+ITSP(J)
ESPAT = ESPAT+ETSP(J)
ENDDO
CSUM=ICSUM
CSUMS=ICSUMS
AVCSUM=CSUM/HN
AVCSMS=CSUMS/HN
AVCDIS=DBLE(ICDI)/HN
CST=DBLE(ICSUM-ICDI)
call writeIntegralSummary(IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,
& ISPAT,ESPAT,tryE,negE,epsilon,prcoeff,L,XX,CW,IRP,
& IRPL,NLI,number_in_layer,depth_interval_flag,AVCSUM,
& AVCDIS,AVCSMS,MAXNL,MAXD)
IF(YH.LE.1.D0) GO TO 1502
AVNLI=ENUCLI/YH
VANLI=ENL2I/YH-AVNLI*AVNLI
SIGNLI=DSQRT(VANLI)
DFINLI=SIGNLI/YH
AVILI=EINELI/YH
VAILI=EIL2I/YH-AVILI*AVILI
SIGILI=DSQRT(VAILI)
DFIILI=SIGILI/YH
CALL MOMENTS(FIX0,SEX,THX,FOX,FIX,SIX,SIGMAX,DFIX0,DSEX,DTHX, XSUM
& ,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,YH)
CALL MOMENTS(FIR0,SER,THR,FOR,FIR,SIR,SIGMAR,DFIR0,DSER,DTHR, RSUM
& ,R2SUM,R3SUM,R4SUM,R5SUM,R6SUM,YH)
CALL MOMENTS(FIP0,SEP,THP,FOP,FIP,SIP,SIGMAP,DFIP0,DSEP,DTHP,
& PLSUM,PL2SUM,PL3SUM,PL4SUM,PL5SUM,PL6SUM,YH)
CALL MOMENTS(FIE0,SEE,THE,FOE,FIE,SIE,SIGMAE,DFIE0,DSEE,DTHE, EEL
& ,EEL2,EEL3,EEL4,EEL5,EEL6,CSUM)
CALL MOMENTS(FIW0,SEW,THW,FOW,FIW,SIW,SIGMAW,DFIW0,DSEW,DTHW,
& EELWC,EELWC2,EELWC3,EELWC4,EELWC5,EELWC6,CSUM)
CALL MOMENTS(FII0,SEI,THI,FOI,FII,SII,SIGMAI,DFII0,DSEI,DTHI, EIL
& ,EIL2,EIL3,EIL4,EIL5,EIL6,CSUM)
CALL MOMENTS(FIS0,SES,THS,FOS,FIS,SIS,SIGMAS,DFIS0,DSES,DTHS, EPL
& ,EPL2,EPL3,EPL4,EPL5,EPL6,CST)
WRITE(21,7117)
7117 FORMAT(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',4X,
& ' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',3X
& ,' ERROR 2.M ', 3X,' ERROR 3.M ')
WRITE(21,7227) FIX0,SEX,THX,FOX,SIGMAX,DFIX0,DSEX,DTHX
7227 FORMAT(1X,' PENETRATION',5X,1P1E12.4,7E14.4)
WRITE(21,7228) FIR0,SER,THR,FOR,SIGMAR,DFIR0,DSER,DTHR
7228 FORMAT(1X,' SPREAD',5X,1P1E12.4,7E14.4)
WRITE(21,7229) FIP0,SEP,THP,FOP,SIGMAP,DFIP0,DSEP,DTHP
7229 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4)
WRITE(21,7237) AVNLI,VANLI,SIGNLI,DFINLI
7237 FORMAT(1X,'ELASTIC LOSS',5X,1P1E12.4,1E14.4,28X,2E14.4)
WRITE(21,7238) AVILI,VAILI,SIGILI,DFIILI
7238 FORMAT(1X,' INEL. LOSS',5X,1P1E12.4,1E14.4,28X,2E14.4)
WRITE(21,7117)
WRITE(21,7231) FIE0,SEE,THE,FOE,SIGMAE,DFIE0,DSEE,DTHE
7231 FORMAT(1X,' ELAST.LOSS',5X,1P1E12.4,7E14.4)
WRITE(21,7232) FIW0,SEW,THW,FOW,SIGMAW,DFIW0,DSEW,DTHW
7232 FORMAT(1X,'WEAK EL.LOSS',5X,1P1E12.4,7E14.4)
WRITE(21,7233) FII0,SEI,THI,FOI,SIGMAI,DFII0,DSEI,DTHI
7233 FORMAT(1X,'INELAST.LOSS',5X,1P1E12.4,7E14.4)
WRITE(21,7234) FIS0,SES,THS,FOS,SIGMAS,DFIS0,DSES,DTHS
7234 FORMAT(1X,' SUBTHR.LOSS',5X,1P1E12.4,7E14.4)
1502 CONTINUE
IF(YH.LT.1.D0) GO TO 7235
CALL MOMENT(X1SD,X2SD,X3SD,X4SD,X5SD,X6SD ,XSUM,X2SUM,X3SUM,X4SUM
& ,X5SUM,X6SUM,YH)
WRITE(21,7118)
WRITE(21,1513) X1SD,X2SD,X3SD,X4SD,X5SD,X6SD
1513 FORMAT(1X,' PENETRATION',5X,1P1E12.4,5E14.4)
7235 continue
if(irl.eq.0) goto 1453
CSUMR=DBLE(ICSUMR)
ACSUMR=CSUMR/HN
ACDISR=DBLE(ICDIR)/HN
ACSBER=DBLE(ICSBR)/HN
WRITE(21,1599) ACSUMR,ACDISR,ACSBER
1599 FORMAT(///2X,'RECOILES (PER PROJ.) : ',
& 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X,
& 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/76X,
& 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3)
IF(NPA+NSA.EQ.0) GO TO 1453
ACSUR=CSUMR/(DBLE(NPA+NSA))
ACDIR=DBLE(ICDIR)/(NPA+NSA)
ACSBR=DBLE(ICSBR)/(NPA+NSA)
WRITE(21,1598) ACSUR,ACDIR,ACSBR
1598 FORMAT(/2X,'RECOILES (PER KNOCKON) : ',
& 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.3,3X,
& 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/,76X,
& 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3/)
IF(NJ(1).LT.2) GO TO 1453
ACDR11=DBLE(ICDIRJ(1,1))/(NPA+NSA)
ACDR12=DBLE(ICDIRJ(1,2))/(NPA+NSA)
ACDR21=DBLE(ICDIRJ(2,1))/(NPA+NSA)
ACDR22=DBLE(ICDIRJ(2,2))/(NPA+NSA)
WRITE(21,1451) ACDR11,ACDR12,ACDR21,ACDR22
1451 FORMAT(76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-1 = ',F10.3/
& ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-2 = ',F10.3/
& ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-1 = ',F10.3/
& ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-2 = ',F10.3)
1453 CONTINUE
call writeImplantationProfile(depth_interval_flag,YH,MAXD,
& MAXD1,LJ,NJ,L,XX,CW,IRL,IRP,RIRP,IPL,ION,DENT,DMGN,
& ELGD,PHON,CASMOT,ICD,ICDT,ICDR,ICDTR,ICDIRI,ICDIRN,
& IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO,
& ICDTT,TIONR,TDENTR,TELGDR,TDMGNR,TPHONR,ICDTTR,ELE,
& ELI,ELP,ELD,ELET,ELIT,ELPT,ELDT,ICDJT,ICDJTR,ICDITR,
& ILD,DLI,IONR,DENTR,DMGNR,ELGDR,PHONR,MAXNL,MAXNL5)
C
C BACKSCATTERING
C
call writeBackscatterSummary(IB,IBL,HN,EB,EMV,ENUCLB,
& ENL2B,EINELB,EIL2B,EB2SUM,EB3SUM,EB4SUM,EB5SUM,
& EB6SUM,EB1SUL,EB2SUL,EB3SUL,EB4SUL,EB5SUL,
& EB6SUL,PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB,
& KADB,AI,RKADB,FIB0,SIGMAB)
C
C TRANSMISSION
C
call writeTransmissionSummary(IT,HN,ET,E0,ENUCLT,ENL2T,
& EINELT,EIL2T,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM,
& PLST,PL2ST,PL3ST,PL4ST,PL5ST,PL6ST,KADT,AI,RKADT,
& FIT0,SIGMAT)
1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3)
C
C BACKWARD SPUTTERING : YIELDS AND ENERGIES
C
call writeBackwardSputterYieldSummary(ISPA,L,NJ,JT,LJ,HN,EMV,
& IBSP,IBSPL,EBSP,ISPAL,ESPAL,ESPA,ISPIP,ISPIS,ISPOP,
& ISPOS,ESPIP,ESPIS,ESPOP,ESPOS,RIP,RIS,ROP,ROS,REIP,
& REIS,REOP,REOS,RIPJ,RISJ,ROPJ,ROSJ,REIPJ,REISJ,
& REOPJ,REOSJ,ESPMIP,ESPMIS,ESPMOP,ESPMOS,SPY,SPE,
& REY,EMSP,SPE2S,SPE3S,SPE4S,SPE5S,SPE6S,SPE1SL,
& SPE2SL,SPE3SL,SPE4SL,SPE5SL,SPE6SL,FIESB,SEESB,
& THESB,FOESB,SGMESB,DFIESB,DSEESB,DTHESB,MAXNL,
& MAXNL5,MAXNLm15)
IF(ISPA.EQ.0) GO TO 1700
1564 FORMAT(9X,'SPUTTERED PARTICLES(',I1,') = ',I7,5X,
& 'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV')
1577 FORMAT(/)
1582 FORMAT(5X,'SPUTTERING YIELD(',I1,') = ',1PE10.3,
& ' SPUTTERED ENERGY(',I1,') = ',1E10.3,
& ' REL.MEAN ENERGY(',I1,') = ',1E10.3,
& ' MEAN ENERGY(',I1,') = ',1E10.3)
7118 FORMAT(/20X,' 1.MOMENT ',4X,' 2.MOMENT ',4X,' 3.MOMENT ',
& 4X,' 4.MOMENT ',4X,' 5.MOMENT ',4X,' 6.MOMENT ')
C
C BACKWARD SPUTTERING : ANGULAR DISTRIBUTIONS
C
call writeBackwardSputterAngularSummary(ISPA,L,NJ,JT,IBSP,
& ISPAL,KADS,KADSJ,KADSL,RKADS,RKADSJ,RKADSL,AI,MAXNL,
& MAXNLm15)
C
C TRANSMISSION SPUTTERING : YIELDS AND ENERGIES
C
1700 IF(ISPAT.EQ.0) GO TO 1800
call writeTransmissionSputterYieldSummary(ISPAT,L,NJ,JT,LJ,NH,
& E0,ESPAT,ITSP,ETSP,ISPALT,ESPALT,ISPIPT,ISPIST,ISPOPT,
& ISPOST,ESPIPT,ESPIST,ESPOPT,ESPOST,RIPT,RIST,ROPT,ROST,
& REIPT,REIST,REOPT,REOST,ESPMIPT,ESPMIST,ESPMOPT,
& ESPMOST,SPYT,SPET,REYT,EMSPT,MAXNL,MAXNL5)
C
C TRANSMISSION SPUTTERING : ANGULAR DISTRIBUTIONS
C
call writeTransmissionSputterAngularSummary(ISPAT,L,NJ,JT,LJ,
& ISPAL,ISPALT,ITSP,KADST,KDSTJ,KDSTL,RKADST,RKDSTJ,
& RKDSTL,AI,MAXNL,MAXNL5,MAXNLm15)
1800 CONTINUE
C
C The per-run summary file (*.seq) is written by writeSequenceSummary.
C
C The *.seq file is the compact per-run summary consumed by the GUI for
C scan plots (fractions stopped in each layer, backscattering,
C transmission, mean depth, etc.). Each run writes its own summary file.
C=======================================================================
C WRITE SEQUENCE SUMMARY FILE
C=======================================================================
call writeSequenceSummary(seqnam,chem,NLayers,E0keV,EsigkeV,
& ALPHA,ALPHASIG,NH,IIM,IB,IT,tryE,negE,FIX0,SIGMAX,FIB0,
& SIGMAB,FIT0,SIGMAT,epsilon,prcoeff,number_in_layer,MAXNL)
C
C End of per-run summary file
C
C TOP AND FRONT LINES FOR MATRICES
C
call setupMatrixTableAxes(DE,DA,DG,NH,NG1,NA1,NGIK,MAXD1,
& JT,MEAB,MAGB,EMA,MEAT,MAGT,EMAT,MEPB,MEPT,MAGS,
& MAGST,MAGSA,MEAS,MEAST,MEASL,MEASTL,MEAGS,MEAGB,
& MEAGT,ELOG,TEMPNH,MAXNL,MAXD2,MAXNLm15)
C
C BACKWARD SPUTTERING : MATRICES , ENERGY - ANGLE CORRELATIONS
C
call writeBackwardSputterMatrixSummary(ISPA,NH,JT,NJ,EASL,
& MEASL,TEMPNH,ELOG,MEAS,MAXD1,ALPHA,NGIK,DGI,MEAGS,
& MAGSA,MAGS,L,ISPAL,MAXNL,MAXD2,MAXNLm15)
C
C Shared matrix output formats used by the remaining inline matrix
C output sections. Copies are kept inside the extracted subroutine.
C
1858 FORMAT(1X,1E12.4,20I5,I6,1E12.4)
1886 FORMAT(1X,I3,20I6,I8)
C
C FORWARD SPUTTERING : MATRICES , ENERGY - ANGLE CORRELATIONS
C
IF(ISPAT.LT.10000) GO TO 2000
JTJ=JT(3)
IF(L.EQ.3) JTJ=LJ-NJ(1)
DO J=1,JTJ
EASTL(2,J)=DBLE(MEASTL(2,21,J))/(DBLE(NH)*0.1D0)
DO IESLOG=3,74
EASTL(IESLOG,J)=DBLE(MEASTL(IESLOG,21,J))/(TEMPNH*10.D0
& **((IESLOG-1)/12.D0))
ENDDO
ENDDO
IF(L.EQ.3) GO TO 1953
DO 1952 J=1,NJ(1)
WRITE(21,1954) J
1954 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(FORWARD SPUTTERED PARTICLES) ',
& ', 1. LAYER , SPECIES',I2/)
do ima = 74,2,-1
if(meastl(ima,21,j).ne.0) goto 1955
enddo
ima = 1
1955 ima = min(ima+2,74)
do ies = 1, ima
write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21),
& eastl(ies,j)
enddo
write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21),eastl(75,j)
WRITE(21,1984) J
1984 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) ',
& '(FORWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/)
do ima = MAXD1,1,-1
if(meast(ima,22,j).ne.0) goto 1983
enddo
ima = 1
1983 ima = min(ima+2,MAXD1)
write (6, 1886) ((meast(iesp,iags,j),iags=1,22),iesp=1,ima)
write (6, 1886) (meast(102,iags,j),iags=1,22)
WRITE(21,1988) J
1988 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) ',
& ' (FORWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/)
WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62)
1952 CONTINUE
1953 CONTINUE
IF(L.EQ.1) GO TO 2000
IF(L.EQ.3) GO TO 1961
JTK=NJ(1)+1
JTL=JT(3)
GO TO 1963
1961 JTK=1
JTL=NJ(2)
1963 DO 1962 J=JTK,JTL
WRITE(21,1964) J-JTK+1
1964 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(FORWARD SPUTTERED PARTICLES) ,',
& ' 2. LAYER , SPECIES',I2/)
do ima = 74,1,-1
if(meastl(ima,21,j).ne.0) goto 1965
enddo
ima = 1
1965 ima = min(ima+2,74)
do ies = 1, ima
write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21),
& eastl(ies,j)
enddo
write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21) ,
& eastl(75,j)
WRITE(21,1994) J-JTK+1
1994 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (FORWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/)
do ima = MAXD1,1,-1
if(meast(ima,22,j).ne.0) goto 1995
enddo
ima = 1
1995 ima = min(ima+2,MAXD1)
WRITE(21,1886)((meast(iesp,iags,j),iags=1,22),iesp=1,ima)
WRITE(21,1886)(meast(102,iags,j),iags=1,22)
WRITE(21,1998) J-JTK+1
1998 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (FORWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/)
WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62)
1962 CONTINUE
IF(L.LT.3) GO TO 2000
DO 1972 J=NJ(2)+1,LJ-NJ(1)
WRITE(21,1974) J-NJ(2)
1974 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(FORWARD SPUTTERED PARTICLES) , 3. LAYER , SPECIES',I2/)
do ima = 74,1,-1
if(meastl(ima,21,j).ne.0) goto 1973
enddo
ima = 1
1973 ima = min(ima+2,74)
do ies = 1, ima
write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21) ,
& eastl(ies,j)
end do
write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21) ,
& eastl(75,j)
WRITE(21,1975) J-NJ(2)
1975 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (FORWARD SPUTTERED PARTICLES) , 3.LAYER , SPECIES',I2/)
do ima = MAXD1,1,-1
if(meast(ima,22,j).ne.0) goto 1977
enddo
ima = 1
1977 ima = min(ima+2,MAXD1)
WRITE(21,1886)((meast(iesp,iags,j),iags=1,22),iesp=1,ima)
WRITE(21,1886)(meast(102,iags,j),iags=1,22)
WRITE(21,1978) J-NJ(2)
1978 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (FORWARD SPUTTERED PARTICLES) , 3.LAYER , SPECIES',I2/)
WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62)
1972 CONTINUE
2000 CONTINUE
C
C BACKSCATTERING : MATRICES , ENERGY - ANGULAR CORRELATIONS
C
IF(IB.LT.10000) GO TO 2100
DO J=1,20
MEABL(1,J)=J
ENDDO
EABL(2)=DBLE(MEABL(2,21))/(DBLE(NH)*0.1D0)
DO IERLOG=3,74
EABL(IERLOG)=DBLE(MEABL(IERLOG,21))/(TEMPNH*10.D0**((IERLOG-1)
& /12.D0))
ENDDO
WRITE(21,2006)
2006 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(BACKSCATTERED PROJECTILES)'/)
do ima = 74,1,-1
if(meabl(ima,21).ne.0) goto 2005
enddo
ima = 1
2005 ima = min(ima+2,74)
do ies = 1, ima
WRITE(21,1858)elog(ies),(meabl(ies,iag),iag=1,21),eabl(ies)
end do
WRITE(21,1858)elog(75),(meabl(75,iag),iag=1,21),eabl(75)
IF(ALPHA.LT.1.) GO TO 2010
DO IG2=1,NGIK,1
EEE = IG2*DGI
WRITE(21,2014) EEE
2014 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) AT AZIMUTHAL ANGLE =',F5.1,
& ' (BACKSCATTERED PROJECTILES)'/)
do ima = MAXD1,1,-1
if(meagb(ima,ig2,22).ne.0) goto 2015
enddo
ima = 1
2015 ima = min(ima+2,MAXD1)
write (6, 1886) ((meagb(ie,ig2,iagb),iagb=1,22),ie=1,ima)
write (6, 1886) (meagb(102,ig2,iagb),iagb=1,22)
ENDDO
2010 CONTINUE
IF(E0.LT.0.) GO TO 2052
WRITE(21,2016)
2016 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (BACKSCATTERED PROJECTILES)'/)
GO TO 2054
2052 WRITE(21,2056)
2056 FORMAT(//,' ENERGY(E IN 0.1*TI) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (BACKSCATTERED PROJECTILES)'/)
do ima = MAXD1,1,-1
if(meab(ima,22).ne.0) goto 2017
enddo
ima = 1
2017 ima = min(ima+2,MAXD1)
write (6, 1886) ((meab(ie,iagb),iagb=1,22),ie=1,ima)
write (6, 1886) (meab(102,iagb),iagb=1,22)
2054 continue
WRITE(21,2018)
2018 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (BACKSCATTERED PROJECTILES)'/)
WRITE(21,1886) ((MAGB(IG,IAGB),IAGB=1,22),IG=1,62)
WRITE(21,2022)
2022 FORMAT(//,' AZIMUTH.ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (BACKSCATTERED ENERGY)'/)
WRITE(21,2027) (EMA(01,IAGB),IAGB=1,11)
WRITE(21,2025) ((EMA(IG,IAGB),IAGB=1,11),IG=2,NG)
WRITE(21,2028)
WRITE(21,2031) EMA(1,1),(EMA(1,IAGB),IAGB=12,22)
DO IG=2,NG
WRITE(21,2026) EMA(IG,1),(EMA(IG,IAGB),IAGB=12,22)
ENDDO
2025 FORMAT(1X,1F5.0,10E11.4)
2026 FORMAT(1X,1F5.0,11E11.4)
2027 FORMAT(1X,1F5.0,10F11.0)
2028 FORMAT(/)
2031 FORMAT(1H1,1X,1F5.0,11F11.0)
2036 FORMAT(1X,26I4)
2038 FORMAT(1X,26I4,I6)
2100 CONTINUE
C
C TRANSMISSION : MATRICES , ENERGY - ANGULAR CORRELATIONS
C
IF(IT.LT.10000) GO TO 9000
DO J=1,20
MEATL(1,J)=J
ENDDO
EATL(2)=DBLE(MEATL(2,21))/(DBLE(NH)*0.1D0)
DO IERLOG=3,74
EATL(IERLOG)=DBLE(MEATL(IERLOG,21))/(TEMPNH* 10.D0**((IERLOG-1)
& /12.D0))
ENDDO
WRITE(21,2106)
2106 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(TRANSMITTED PROJECTILES)'/)
do ima = 74,1,-1
if(meatl(ima,21).ne.0) goto 2105
enddo
ima = 1
2105 ima = min(ima+2,74)
do ies = 1, ima
WRITE(21,1858)elog(ies),(meatl(ies,iag),iag=1,21),eatl(ies)
enddo
WRITE(21,1858)elog(75),(meatl(75,iag),iag=1,21),eatl(75)
IF(ALPHA.LT.1.) GO TO 2110
DO IG2=1,NGIK,1
EEE = IG2*DGI
WRITE(21,2114) EEE
2114 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) AT AZIMUTHAL ANGLE =',F5.1,
& ' (TRANSMITTED PROJECTILES)'/)
do ima = MAXD1,1,-1
if(meagt(ima,ig2,22).ne.0) goto 2115
enddo
ima = 1
2115 ima = min(ima+2,MAXD1)
write (21,1886) ((meagt(ie,ig2,iagb),iagb=1,22),ie=1,ima)
write (21,1886) (meagt(102,ig2,iagb),iagb=1,22)
ENDDO
2110 CONTINUE
WRITE(21,2116)
2116 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (TRANSMITTED PROJECTILES)'/)
do ima = MAXD1,1,-1
if(meat(ima,22).ne.0) goto 2117
enddo
ima = 1
2117 ima = min(ima+2,MAXD1)
write (6, 1886) ((meat(ie,iagb),iagb=1,22),ie=1,ima)
write (6, 1886) (meat(102,iagb),iagb=1,22)
WRITE(21,2118)
2118 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (TRANSMITTED PROJECTILES)'/)
WRITE(21,1886) ((MAGT(IG,IAGB),IAGB=1,22),IG=1,62)
WRITE(21,2122)
2122 FORMAT(//,' AZIMUTH.ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (TRANSMITTED ENERGY)'/)
WRITE(21,2127) (EMAT(01,IAGB),IAGB=1,11)
WRITE(21,2125) ((EMAT(IG,IAGB),IAGB=1,11),IG=2,NG)
WRITE(21,2028)
WRITE(21,2131) EMAT(1,1),(EMAT(1,IAGB),IAGB=12,22)
DO IG=2,NG
WRITE(21,2126) EMAT(IG,1),(EMAT(IG,IAGB),IAGB=12,22)
ENDDO
2125 FORMAT(1X,1F5.0,10E11.4)
2126 FORMAT(1X,1F5.0,11E11.4)
2127 FORMAT(1X,1F5.0,10F11.0)
2131 FORMAT(1H1,1X,1F5.0,11F11.0)
9000 CONTINUE
CLOSE(UNIT=21)
CLOSE(UNIT=22)
CLOSE(UNIT=99)
8000 CONTINUE
END
C=======================================================================
C OUTPUT WRITER SUBROUTINES
C=======================================================================
C
C Write backward sputtering matrix and energy-angle correlation tables.
C
subroutine writeBackwardSputterMatrixSummary(ISPA,NH,JT,NJ,
& EASL,MEASL,TEMPNH,ELOG,MEAS,MAXD1,ALPHA,NGIK,DGI,
& MEAGS,MAGSA,MAGS,L,ISPAL,MAXNL,MAXD2,MAXNLm15)
implicit none
integer*4 ISPA,NH,MAXD1,NGIK,L,MAXNL,MAXD2,MAXNLm15
integer*4 J,IESLOG,IMA,IES,IAS,IESP,IAGS,IG2,IG,IA
integer*4 JT(MAXNL),NJ(MAXNL),ISPAL(MAXNL)
integer*4 MEASL(75,21,MAXNLm15),MEAS(MAXD2,22,MAXNLm15)
integer*4 MEAGS(MAXD2,12,22,MAXNLm15)
integer*4 MAGSA(62,32,MAXNLm15),MAGS(62,22,MAXNLm15)
real*8 EASL(75,MAXNLm15),TEMPNH,ELOG(75),ALPHA,DGI,EEE
IF(ISPA.LT.10000) GO TO 1900
DO J=1,JT(3)
EASL(2,J)=DBLE(MEASL(2,21,J))/(DBLE(NH)*0.1)
DO IESLOG=3,74
EASL(IESLOG,J)=DBLE(MEASL(IESLOG,21,J))/(TEMPNH*10.D0
& **((IESLOG-1)/12.D0))
ENDDO
ENDDO
DO 1852 J=1,NJ(1)
WRITE(21,1854) J
1854 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (BAC
&KWARD SPUTTERED PARTICLES) ; 1. LAYER ; SPECIES',I2/)
do ima = 74,2,-1
if(measl(ima,21,j).ne.0) goto 1855
enddo
ima = 1
1855 ima = min(ima+2,74)
do ies = 1, ima
write (6, 1858) elog(ies), (measl(ies,ias,j),ias=1,21),
& easl(ies,j)
enddo
write (6, 1858) elog(75), (measl(75,ias,j),ias=1,21),easl(75,j)
1858 FORMAT(1X,1E12.4,20I5,I6,1E12.4)
WRITE(21,1884) J
1884 FORMAT(//,' ENERGY(E/E0 IN %) - ',
& 'POLAR ANGLE IN COS-INTERVALS (0.05) ',
& '(BACKWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/)
do ima = MAXD1,1,-1
if(meas(ima,22,j).ne.0) goto 1883
enddo
ima = 1
1883 ima = min(ima+2,MAXD1)
write (6, 1886) ((meas(iesp,iags,j),iags=1,22),iesp=1,ima)
write (6, 1886) (meas(102,iags,j),iags=1,22)
1886 FORMAT(1X,I3,20I6,I8)
IF(ALPHA.LT.1.) GO TO 1878
DO 1870 IG2=1,NGIK,1
EEE = IG2*DGI
WRITE(21,1872) EEE,J
1872 FORMAT(//,' ENERGY(E/E0 IN %) - ',
& 'POLAR ANGLE IN COS-INTERVALS (0.05) ',
& 'AT AZIMUTHAL ANGLE =',F5.1,
& ' (BACKWARD SPUTTERED ATOMS) , 1.LAYER , SPECIES',I2/)
do ima = MAXD1,1,-1
if(meags(ima,ig2,22,j).ne.0) goto 1885
enddo
ima = 1
1885 ima = min(ima+2,MAXD1)
do iesp = 1, ima
write (6, 1886) (meags(iesp,ig2,iags,j),iags=1,22)
end do
write (6, 1886) (meags(102,ig2,iags,j),iags=1,22)
1870 CONTINUE
WRITE(21,1889) J
1889 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN DEGREES ',
& '(BACKWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/)
WRITE(21,1887) ((MAGSA(IG,IA,J),IA=1,32),IG=1,62)
1887 FORMAT(1X,31I4,I6)
1878 CONTINUE
WRITE(21,1888) J
1888 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) ',
& ' (BACKWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/)
WRITE(21,1886) ((MAGS(IG,IAGS,J),IAGS=1,22),IG=1,62)
1852 CONTINUE
IF(L.EQ.1) GO TO 1900
if(ispal(2).eq.0) goto 1900
DO 1862 J=NJ(1)+1,JT(3)
WRITE(21,1864) J-NJ(1)
1864 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) ',
& '(BACKWARD SPUTTERED PARTICLES) , 2. LAYER , SPECIES',I2/)
do ima = 74,1,-1
if(measl(ima,21,j).ne.0) goto 1865
enddo
ima = 1
1865 ima = min(ima+2,74)
do ies = 1, ima
write (6, 1858) elog(ies),(measl(ies,ias,j),ias=1,21),
& easl(ies,j)
enddo
write (6, 1858) elog(75),(measl(75,ias,j),ias=1,21),easl(75,j)
WRITE(21,1894) J-NJ(1)
1894 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) (BACKWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',
& I2/)
do ima = MAXD1,1,-1
if(meas(ima,22,j).ne.0) goto 1895
enddo
ima = 1
1895 ima = min(ima+2,MAXD1)
WRITE(21,1886)((meas(iesp,iags,j),iags=1,22),iesp=1,ima)
WRITE(21,1886)(meas(102,iags,j),iags=1,22)
WRITE(21,1898) J-NJ(1)
1898 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS ',
& '(0.05) ',
1 ' (BACKWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/)
WRITE(21,1886) ((MAGS(IG,IAGS,J),IAGS=1,22),IG=1,62)
1862 CONTINUE
1900 CONTINUE
return
end
C
C Write angular distributions for transmission sputtered particles.
C
subroutine writeTransmissionSputterAngularSummary(ispat,nLayers,
& nj,jt,nSpecies,ispal,ispalt,itsp,kadst,kdstj,kdstl,
& rkadst,rkdstj,rkdstl,angleBins,maxNl,maxNl5,maxNlm15)
implicit none
integer maxNl,maxNl5,maxNlm15
integer ispat,nLayers,nSpecies
integer i,j
integer*4 nj(maxNl),jt(maxNl)
integer*4 ispal(maxNl),ispalt(maxNl),itsp(maxNl5)
integer*4 kadst(20),kdstj(20,maxNlm15),kdstl(20,6)
real*8 rkadst(20),rkdstj(20,maxNlm15),rkdstl(20,maxNl)
real*8 angleBins(20)
write(21,1760)
1760 format(///5X,'POLAR ANGULAR DISTRIBUTION OF ALL TRANSMISSION '//
& 'SPUTTERED PARTICLES'//)
do i=1,20
rkadst(i)=kadst(i)*20.D0/ispat
enddo
write(21,1518) (angleBins(i),i=1,20),(kadst(i),i=1,20),
& (rkadst(i),i=1,20)
if(nLayers.EQ.3) go to 1764
do i=1,20
do j=1,nj(1)
kdstl(i,1)=kdstl(i,1)+kdstj(i,j)
enddo
enddo
do i=1,20
do j=nj(1)+1,jt(3)
kdstl(i,2)=kdstl(i,2)+kdstj(i,j)
enddo
enddo
1766 continue
do j=1,2
if(ispal(j).EQ.0) go to 1754
do i=1,20
rkdstl(i,j)=kdstl(i,j)*20.D0/ispal(j)
enddo
1754 continue
enddo
do j=1,jt(3)
if(itsp(j).EQ.0) go to 1756
do i=1,20
rkdstj(i,j)=kdstj(i,j)*20.D0/itsp(j)
enddo
1756 continue
enddo
write(21,1776)
1776 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 1'//)
write(21,1518) (angleBins(i),i=1,20),(kdstl(i,1),i=1,20)
& ,(rkdstl(i,1),i=1,20)
if(nj(1).EQ.1) go to 1778
do 1780 j=1,nj(1)
if(itsp(j).EQ.0) go to 1780
write(21,1782) j
1782 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 1 , SPECIES ',I1//)
write(21,1518) (angleBins(i),i=1,20),(kdstj(i,j),i=1,20)
& ,(rkdstj(i,j),i=1,20)
1780 continue
1778 if(nLayers.EQ.1) go to 1800
write(21,1786)
1786 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2'//)
write(21,1518) (angleBins(i),i=1,20),(kdstl(i,2),i=1,20)
& ,(rkdstl(i,2),i=1,20)
if(nj(2).EQ.1) go to 1800
do j=nj(1)+1,jt(3)
write(21,1790) j-nj(1)
1790 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2 , SPECIES ',I1//)
write(21,1518) (angleBins(i),i=1,20),(kdstj(i,j),i=1,20)
& ,(rkdstj(i,j),i=1,20)
enddo
go to 1800
1764 do i=1,20
do j=1,nj(2)
kdstl(i,1)=kdstl(i,1)+kdstj(i,j)
enddo
do j=nj(2)+1,nSpecies-nj(1)
kdstl(i,2)=kdstl(i,2)+kdstj(i,j)
enddo
enddo
do 1799 j=1,2
if(ispalt(j+1).EQ.0) go to 1799
do i=1,20
rkdstl(i,j)=kdstl(i,j)*20.D0/ispalt(j+1)
enddo
1799 continue
do 1797 j=1,nSpecies-nj(1)
if(itsp(j+nj(1)).EQ.0) go to 1797
do i=1,20
rkdstj(i,j)=kdstj(i,j)*20.D0/itsp(j+nj(1))
enddo
1797 continue
write(21,1771)
1771 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2'//)
write(21,1518) (angleBins(i),i=1,20),(kdstl(i,1),i=1,20)
& ,(rkdstl(i,1),i=1,20)
if(nj(2).EQ.1) go to 1773
do 1775 j=1,nj(2)
if(itsp(j+nj(1)).EQ.0) go to 1775
write(21,1777) j
1777 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2 ; SPECIES ',I1//)
write(21,1518) (angleBins(i),i=1,20),(kdstj(i,j),i=1,20)
& ,(rkdstj(i,j),i=1,20)
1775 continue
1773 continue
write(21,1779)
1779 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 3'//)
write(21,1518) (angleBins(i),i=1,20),(kdstl(i,2),i=1,20)
& ,(rkdstl(i,2),i=1,20)
if(nj(2).EQ.1) go to 1800
do j=nj(2)+1,nSpecies-nj(1)
write(21,1783) j-nj(2)
1783 format(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 3 ; SPECIES ',I1//)
write(21,1518) (angleBins(i),i=1,20),(kdstj(i,j),i=1,20)
& ,(rkdstj(i,j),i=1,20)
enddo
1800 continue
return
1518 format(5X,20F6.2//,5X,20I6/5X,20F6.3)
end
C
C Write transmission sputtering yields and energy summaries.
C
subroutine writeTransmissionSputterYieldSummary(ispat,nLayers,
& nj,jt,nSpecies,nProjectiles,initialEnergy,totalEnergy,
& itsp,etsp,ispalt,espalt,ispipt,ispist,ispopt,ispost,
& espipt,espist,espopt,espost,ript,rist,ropt,rost,reipt,
& reist,reopt,reost,espmipt,espmist,espmopt,espmost,spyt,
& spet,reyt,emspt,maxNl,maxNl5)
implicit none
integer maxNl,maxNl5
integer ispat,nLayers,nSpecies,nProjectiles
integer j
integer*4 nj(maxNl),jt(maxNl)
integer*4 itsp(maxNl5),ispalt(maxNl)
integer*4 ispipt(maxNl5),ispist(maxNl5)
integer*4 ispopt(maxNl5),ispost(maxNl5)
real*8 initialEnergy,totalEnergy
real*8 etsp(maxNl5),espalt(maxNl)
real*8 ript(maxNl5),rist(maxNl5),ropt(maxNl5),rost(maxNl5)
real*8 reipt(maxNl5),reist(maxNl5),reopt(maxNl5)
real*8 reost(maxNl5)
real*8 espipt(maxNl5),espist(maxNl5),espopt(maxNl5)
real*8 espost(maxNl5)
real*8 espmipt(maxNl5),espmist(maxNl5),espmopt(maxNl5)
real*8 espmost(maxNl5)
real*8 spyt(maxNl5),spet(maxNl5),reyt(maxNl5),emspt(maxNl5)
logical equal
external equal
if(ispat.eq.0) return
write(21,1704)
1704 format(1H1,5X,'TRANSMISSION SPUTTERING')
do j=1,nj(1)
ispalt(1) = ispalt(1)+itsp(j)
espalt(1) = espalt(1)+etsp(j)
enddo
do j=nj(1)+1,jt(3)
ispalt(2) = ispalt(2)+itsp(j)
espalt(2) = espalt(2)+etsp(j)
enddo
do j=jt(3)+1,nSpecies
ispalt(3) = ispalt(3)+itsp(j)
espalt(3) = espalt(3)+etsp(j)
enddo
write(21,1712) ispat,totalEnergy
1712 format(///,8X,'ALL SPUTTERED PARTICLES = ',I7,3X,
& 'TOTAL SPUTTERED ENERGY = ',E10.4,3H EV//)
do j=1,nLayers
write(21,1713) j,ispalt(j),espalt(j)
1713 format(8X,'SPUTTERED PARTICLES (LAYER ',I3,') = ',I7,
& 3X,'SPUTTERED ENERGY = ',E10.4,3H EV)
enddo
do j=1,nSpecies
ript(j)=dble(ispipt(j))/dble(ispat)
rist(j)=dble(ispist(j))/dble(ispat)
ropt(j)=dble(ispopt(j))/dble(ispat)
rost(j)=dble(ispost(j))/dble(ispat)
reipt(j)=espipt(j)/totalEnergy
reist(j)=espist(j)/totalEnergy
reopt(j)=espopt(j)/totalEnergy
reost(j)=espost(j)/totalEnergy
enddo
1715 continue
do 1717 j=1,nSpecies
if(ispipt(j).eq.0) go to 4571
espmipt(j)=espipt(j)/dble(ispipt(j))
4571 if(ispist(j).eq.0) go to 4572
espmist(j)=espist(j)/dble(ispist(j))
4572 if(ispopt(j).eq.0) go to 4573
espmopt(j)=espopt(j)/dble(ispopt(j))
4573 if(ispost(j).eq.0) go to 1717
espmost(j)=espost(j)/dble(ispost(j))
1717 continue
do j=1,nSpecies
spyt(j)=dble(itsp(j))/dble(nProjectiles)
spet(j)=etsp(j)/(nProjectiles*initialEnergy)
enddo
do 1737 j=1,nSpecies
if (equal(spyt(j),0.0D0))go to 1737
reyt(j)=spet(j)/spyt(j)
emspt(j)=reyt(j)*initialEnergy
1737 continue
if(ispalt(1).eq.0) go to 1719
write(21,1714)
1714 format(//1X,'1.LAYER')
do j=1,nj(1)
write(21,1564) j,itsp(j),j,etsp(j)
enddo
do 1734 j=1,nj(1)
write(21,1581) j,ispipt(j),ript(j),espipt(j),reipt(j),
& espmipt(j),j,ispist(j),rist(j),espist(j),reist(j),
& espmist(j),j,ispopt(j),ropt(j),espopt(j),reopt(j),
& espmopt(j),j,ispost(j),rost(j),espost(j),reost(j),
& espmost(j)
1734 continue
1581 format(/9X,'ION IN , PRIMARY KO(',I1,') = ',I7,1F9.4,
& 4X,'ENERGY = ',E10.4,' EV',1F9.4,4X,
& 'MEAN ENERGY = ',E10.4/
& 9X,'ION IN , SECOND. KO(',I1,') = ',I7,1F9.4,
& 4X,'ENERGY = ',E10.4,' EV',1F9.4,4X,
& 'MEAN ENERGY = ',E10.4/
& 8X,'ION OUT , PRIMARY KO(',I1,') = ',I7,1F9.4,
& 4X,'ENERGY = ',E10.4,' EV',1F9.4,4X,
& 'MEAN ENERGY = ',E10.4/
& 8X,'ION OUT , SECOND. KO(',I1,') = ',I7,1F9.4,
& 4X,'ENERGY = ',E10.4,' EV',1F9.4,4X,
& 'MEAN ENERGY = ',E10.4)
write(21,1577)
do j=1,nj(1)
write(21,1582) j,spyt(j),j,spet(j),j,reyt(j),j,emspt(j)
enddo
1719 if(nLayers.eq.1) go to 1749
if(ispalt(2).eq.0) go to 1744
write(21,1720)
1720 format(/1X,'2.LAYER')
do j=nj(1)+1,jt(3)
write(21,1564) j-nj(1),itsp(j),j-nj(1),etsp(j)
enddo
do j=nj(1)+1,jt(3)
write(21,1581) j-nj(1),ispipt(j),ript(j),espipt(j),
& reipt(j),espmipt(j),j-nj(1),ispist(j),rist(j),
& espist(j),reist(j),espmist(j),j-nj(1),ispopt(j),
& ropt(j),espopt(j),reopt(j),espmopt(j),j-nj(1),
& ispost(j),rost(j),espost(j),reost(j),espmost(j)
enddo
write(21,1577)
do j=nj(1)+1,jt(3)
write(21,1582) j-nj(1),spyt(j),j-nj(1),spet(j),
& j-nj(1),reyt(j),j-nj(1),emspt(j)
enddo
1744 if(nLayers.eq.2) go to 1749
if(ispalt(3).eq.0) go to 1749
write(21,1726)
1726 format(/1X,'3.LAYER')
do j=jt(3)+1,nSpecies
write(21,1564) j-jt(3),itsp(j),j-jt(3),etsp(j)
enddo
do j=jt(3)+1,nSpecies
write(21,1581) j-jt(3),ispipt(j),ript(j),espipt(j),
& reipt(j),espmipt(j),j-jt(3),ispist(j),rist(j),
& espist(j),reist(j),espmist(j),j-jt(3),ispopt(j),
& ropt(j),espopt(j),reopt(j),espmopt(j),j-jt(3),
& ispost(j),rost(j),espost(j),reost(j),espmost(j)
enddo
write(21,1577)
do j=jt(3)+1,nSpecies
write(21,1582) j-jt(3),spyt(j),j-jt(3),spet(j),
& j-jt(3),reyt(j),j-jt(3),emspt(j)
enddo
1749 continue
return
1564 format(9X,'SPUTTERED PARTICLES(',I1,') = ',I7,5X,
& 'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV')
1577 format(/)
1582 format(5X,'SPUTTERING YIELD(',I1,') = ',1PE10.3,
& ' SPUTTERED ENERGY(',I1,') = ',1E10.3,
& ' REL.MEAN ENERGY(',I1,') = ',1E10.3,
& ' MEAN ENERGY(',I1,') = ',1E10.3)
end
C=======================================================================
C=======================================================================
C WRITE BACKWARD SPUTTERING ANGULAR SUMMARY
C=======================================================================
subroutine writeBackwardSputterAngularSummary(ISPA,L,NJ,JT,IBSP,
& ISPAL,KADS,KADSJ,KADSL,RKADS,RKADSJ,RKADSL,AI,MAXNL,
& MAXNLm15)
implicit none
integer ISPA,L,MAXNL,MAXNLm15
integer NJ(MAXNL),JT(MAXNL)
integer IBSP(*),ISPAL(MAXNL)
integer KADS(20),KADSJ(20,MAXNLm15),KADSL(20,6)
integer I,J
real*8 AI(20),RKADS(20),RKADSJ(20,MAXNLm15)
real*8 RKADSL(20,MAXNL)
C
WRITE(21,1601)
1601 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF ALL BACKWARD ',
& 'SPUTTERED PARTICLES'//)
DO I=1,20
RKADS(I)=KADS(I)*20.D0/ISPA
ENDDO
WRITE(21,1518) (AI(I),I=1,20),(KADS(I),I=1,20),
& (RKADS(I),I=1,20)
DO I=1,20
DO J=1,NJ(1)
KADSL(I,1)=KADSL(I,1)+KADSJ(I,J)
ENDDO
ENDDO
DO I=1,20
DO J=NJ(1)+1,JT(3)
KADSL(I,2)=KADSL(I,2)+KADSJ(I,J)
ENDDO
ENDDO
IF(ISPAL(1).EQ.0) GO TO 1614
IF(NJ(1).EQ.1) GO TO 1614
WRITE(21,1606)
1606 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 1'//)
DO I=1,20
RKADSL(I,1)=KADSL(I,1)*20.D0/ISPAL(1)
ENDDO
WRITE(21,1518) (AI(I),I=1,20),(KADSL(I,1),I=1,20),
& (RKADSL(I,1),I=1,20)
DO 1618 J=1,NJ(1)
IF(IBSP(J).EQ.0) GO TO 1618
WRITE(21,1616) J
1616 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 1 ; SPECIES ',I1//)
DO I=1,20
RKADSJ(I,J)=KADSJ(I,J)*20.D0/IBSP(J)
ENDDO
WRITE(21,1518) (AI(I),I=1,20),(KADSJ(I,J),I=1,20),
& (RKADSJ(I,J),I=1,20)
1618 CONTINUE
1614 IF(L.EQ.1) GO TO 1622
IF(ISPAL(2).EQ.0) GO TO 1622
WRITE(21,1610)
1610 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2'//)
DO I=1,20
RKADSL(I,2)=KADSL(I,2)*20.D0/ISPAL(2)
ENDDO
WRITE(21,1518) (AI(I),I=1,20),(KADSL(I,2),I=1,20),
& (RKADSL(I,2),I=1,20)
IF(NJ(2).EQ.1) GO TO 1622
DO 1624 J=NJ(1)+1,JT(3)
IF(IBSP(J).EQ.0) GO TO 1624
WRITE(21,1626) J-NJ(1)
1626 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ',
& 'PARTICLES ; LAYER 2 ; SPECIES ',I1//)
DO I=1,20
RKADSJ(I,J)=KADSJ(I,J)*20.D0/IBSP(J)
ENDDO
WRITE(21,1518) (AI(I),I=1,20),(KADSJ(I,J),I=1,20),
& (RKADSJ(I,J),I=1,20)
1624 CONTINUE
1622 CONTINUE
1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3)
return
end
C=======================================================================
C WRITE BACKWARD SPUTTERING YIELD SUMMARY
C=======================================================================
subroutine writeBackwardSputterYieldSummary(ISPA,L,NJ,JT,LJ,HN,
& EMV,IBSP,IBSPL,EBSP,ISPAL,ESPAL,ESPA,ISPIP,ISPIS,ISPOP,
& ISPOS,ESPIP,ESPIS,ESPOP,ESPOS,RIP,RIS,ROP,ROS,REIP,
& REIS,REOP,REOS,RIPJ,RISJ,ROPJ,ROSJ,REIPJ,REISJ,
& REOPJ,REOSJ,ESPMIP,ESPMIS,ESPMOP,ESPMOS,SPY,SPE,
& REY,EMSP,SPE2S,SPE3S,SPE4S,SPE5S,SPE6S,SPE1SL,
& SPE2SL,SPE3SL,SPE4SL,SPE5SL,SPE6SL,FIESB,SEESB,
& THESB,FOESB,SGMESB,DFIESB,DSEESB,DTHESB,MAXNL,
& MAXNL5,MAXNLm15)
implicit none
integer ISPA,L,LJ,MAXNL,MAXNL5,MAXNLm15
integer NJ(MAXNL),JT(MAXNL)
integer IBSP(MAXNL5),IBSPL(MAXNL5),ISPAL(MAXNL)
integer ISPIP(MAXNL5),ISPIS(MAXNL5),ISPOP(MAXNL5)
integer ISPOS(MAXNL5)
integer J
real*8 HN,EMV,ESPA
real*8 EBSP(MAXNL5),ESPAL(MAXNL)
real*8 ESPIP(MAXNL5),ESPIS(MAXNL5),ESPOP(MAXNL5)
real*8 ESPOS(MAXNL5)
real*8 RIP(MAXNL5),RIS(MAXNL5),ROP(MAXNL5),ROS(MAXNL5)
real*8 REIP(MAXNL5),REIS(MAXNL5),REOP(MAXNL5),REOS(MAXNL5)
real*8 RIPJ(MAXNL5),RISJ(MAXNL5),ROPJ(MAXNL5),ROSJ(MAXNL5)
real*8 REIPJ(MAXNL5),REISJ(MAXNL5),REOPJ(MAXNL5)
real*8 REOSJ(MAXNL5)
real*8 ESPMIP(MAXNL5),ESPMIS(MAXNL5),ESPMOP(MAXNL5)
real*8 ESPMOS(MAXNL5)
real*8 SPY(MAXNL5),SPE(MAXNL5),REY(MAXNL5),EMSP(MAXNL5)
real*8 SPE2S(MAXNL5),SPE3S(MAXNL5),SPE4S(MAXNL5)
real*8 SPE5S(MAXNL5),SPE6S(MAXNL5)
real*8 SPE1SL(MAXNL5),SPE2SL(MAXNL5),SPE3SL(MAXNL5)
real*8 SPE4SL(MAXNL5),SPE5SL(MAXNL5),SPE6SL(MAXNL5)
real*8 FIESB(MAXNLm15),SEESB(MAXNLm15),THESB(MAXNLm15)
real*8 FOESB(MAXNLm15),SGMESB(MAXNLm15),DFIESB(MAXNLm15)
real*8 DSEESB(MAXNLm15),DTHESB(MAXNLm15)
real*8 FIES0,SEES,THES,FOES,FIES,SIES,SIGMES
real*8 DFIES0,DSEES,DTHES
real*8 FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,SIGMSL
real*8 DFIESL,DSEESL,DTHESL
real*8 EBSP1,EBSP2,EBSP3,EBSP4,EBSP5,EBSP6
real*8 EBSP1L,EBSP2L,EBSP3L,EBSP4L,EBSP5L,EBSP6L
real*8 YSP,YSPL
logical equal
if(ISPA.eq.0) return
write(21,1548)
1548 format(1H1,5X,'BACKWARD SPUTTERING')
do J=1,NJ(1)
ISPAL(1) = ISPAL(1)+IBSP(J)
ESPAL(1) = ESPAL(1)+EBSP(J)
enddo
do J=NJ(1)+1,JT(3)
ISPAL(2) = ISPAL(2)+IBSP(J)
ESPAL(2) = ESPAL(2)+EBSP(J)
enddo
do J=JT(3)+1,LJ
ISPAL(3) = ISPAL(3)+IBSP(J)
ESPAL(3) = ESPAL(3)+EBSP(J)
enddo
write(21,1558) ISPA,ESPA
1558 format(///,8X,'ALL SPUTTERED PARTICLES = ',I7,3X,
& 'TOTAL SPUTTERED ENERGY = ',E10.4,3H EV//)
do 1557 J=1,L
write(21,1559) J,ISPAL(J),ESPAL(J)
1559 format(8X,'SPUTTERED PARTICLES (',I3,'.LAYER) = ',I7,3X,
& 'SPUTTERED ENERGY = ',E10.4,3H EV)
1557 continue
write(21,1560)
1560 format(//1X,'1.LAYER')
do 1562 J=1,NJ(1)
write(21,1564) J,IBSP(J),J,EBSP(J)
1564 format(9X,'SPUTTERED PARTICLES(',I1,') = ',I7,5X,
& 'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV')
1562 continue
if(ISPA.eq.0) return
do J=1,LJ
RIP(J)=dble(ISPIP(J))/dble(ISPA)
RIS(J)=dble(ISPIS(J))/dble(ISPA)
ROP(J)=dble(ISPOP(J))/dble(ISPA)
ROS(J)=dble(ISPOS(J))/dble(ISPA)
REIP(J)=ESPIP(J)/ESPA
REIS(J)=ESPIS(J)/ESPA
REOP(J)=ESPOP(J)/ESPA
REOS(J)=ESPOS(J)/ESPA
enddo
do 1584 J=1,LJ
if(IBSP(J).eq.0) go to 1584
RIPJ(J)=dble(ISPIP(J))/dble(IBSP(J))
RISJ(J)=dble(ISPIS(J))/dble(IBSP(J))
ROPJ(J)=dble(ISPOP(J))/dble(IBSP(J))
ROSJ(J)=dble(ISPOS(J))/dble(IBSP(J))
REIPJ(J)=ESPIP(J)/EBSP(J)
REISJ(J)=ESPIS(J)/EBSP(J)
REOPJ(J)=ESPOP(J)/EBSP(J)
REOSJ(J)=ESPOS(J)/EBSP(J)
1584 continue
do 1571 J=1,LJ
if(ISPIP(J).eq.0) go to 3571
ESPMIP(J)=ESPIP(J)/dble(ISPIP(J))
3571 if(ISPIS(J).eq.0) go to 3572
ESPMIS(J)=ESPIS(J)/dble(ISPIS(J))
3572 if(ISPOP(J).eq.0) go to 3573
ESPMOP(J)=ESPOP(J)/dble(ISPOP(J))
3573 if(ISPOS(J).eq.0) go to 1571
ESPMOS(J)=ESPOS(J)/dble(ISPOS(J))
1571 continue
1573 continue
do J=1,LJ
SPY(J)=dble(IBSP(J))/HN
SPE(J)=EBSP(J)/(HN*EMV)
enddo
do 1579 J=1,LJ
if (equal(SPY(J),0.0D0)) go to 1579
REY(J)=SPE(J)/SPY(J)
EMSP(J)=EBSP(J)/IBSP(J)
1579 continue
if(ISPAL(1).eq.0) go to 1575
do J=1,NJ(1)
write(21,1576) J,ISPIP(J),RIP(J),RIPJ(J),ESPIP(J),REIP(J),
& REIPJ(J),ESPMIP(J),J,ISPIS(J),RIS(J),RISJ(J),ESPIS(J),
& REIS(J),REISJ(J),ESPMIS(J),J,ISPOP(J),ROP(J),ROPJ(J),
& ESPOP(J),REOP(J),REOPJ(J),ESPMOP(J),J,ISPOS(J),ROS(J),
& ROSJ(J),ESPOS(J),REOS(J),REOSJ(J),ESPMOS(J)
1576 format(/9X,'ION IN , PRIMARY KO(',I1,') = ',I7,2F9.4,4X,
& 'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/
& 9X,'ION IN , SECOND. KO(',I1,') = ',I7,2F9.4,4X,
& 'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/
& 8X,'ION OUT , PRIMARY KO(',I1,') = ',I7,2F9.4,4X,
& 'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/
& 8X,'ION OUT , SECOND. KO(',I1,') = ',I7,2F9.4,4X,
& 'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4)
enddo
1575 continue
write(21,1577)
1577 format(/)
do 1580 J=1,NJ(1)
write(21,1582) J,SPY(J),J,SPE(J),J,REY(J),J,EMSP(J)
1582 format(5X,'SPUTTERING YIELD(',I1,') = ',1PE10.3,
& ' SPUTTERED ENERGY(',I1,') = ',1E10.3,
& ' REL.MEAN ENERGY(',I1,') = ',1E10.3,
& ' MEAN ENERGY(',I1,') = ',1E10.3)
1580 continue
do 7260 J=1,NJ(1)
if(IBSP(J).le.1) go to 7260
YSP=IBSP(J)
YSPL=IBSPL(J)
call MOMENTN(FIES0,SEES,THES,FOES,FIES,SIES,SIGMES,
& DFIES0,DSEES,DTHES,EBSP1,EBSP2,EBSP3,EBSP4,EBSP5,
& EBSP6,EBSP(J),SPE2S(J),SPE3S(J),SPE4S(J),SPE5S(J),
& SPE6S(J),YSP)
call MOMENTN(FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,
& SIGMSL,DFIESL,DSEESL,DTHESL,EBSP1L,EBSP2L,EBSP3L,
& EBSP4L,EBSP5L,EBSP6L,SPE1SL(J),SPE2SL(J),SPE3SL(J),
& SPE4SL(J),SPE5SL(J),SPE6SL(J),YSPL)
write(21,7117)
write(21,7261) J,FIES0,SEES,THES,FOES,SIGMES,DFIES0,
& DSEES,DTHES
7261 format(1X,' ENERGY(',I1,')',5X,1P1E12.4,7E14.4)
write(21,7263) J,FIES0L,SEESL,THESL,FOESL,SIGMSL,
& DFIESL,DSEESL,DTHESL
7263 format(1X,'LOGENERGY(',I1,')',5X,1P1E12.4,7E14.4)
write(21,7118)
7118 format(/20X,' 1.MOMENT ',4X,' 2.MOMENT ',4X,' 3.MOMENT ',
& 4X,' 4.MOMENT ',4X,' 5.MOMENT ',4X,' 6.MOMENT ')
write(21,7265) J,EBSP1,EBSP2,EBSP3,EBSP4,EBSP5,EBSP6
7265 format(1X,' ENERGY(',I1,')',5X,1P1E12.4,5E14.4)
write(21,7267) J,EBSP1L,EBSP2L,EBSP3L,EBSP4L,EBSP5L,
& EBSP6L
7267 format(1X,'LOGENERGY(',I1,')',5X,1P1E12.4,5E14.4)
FIESB(J)=FIES0
SEESB(J)=SEES
THESB(J)=THES
FOESB(J)=FOES
SGMESB(J)=SIGMES
DFIESB(J)=DFIES0
DSEESB(J)=DSEES
DTHESB(J)=DTHES
7260 continue
if(L.eq.1) go to 1593
if(ISPAL(2).eq.0) go to 1593
write(21,1566)
1566 format(//1X,'2.LAYER')
do 1568 J=NJ(1)+1,JT(3)
write(21,1570) J-NJ(1),IBSP(J),J-NJ(1),EBSP(J)
1570 format(9X,'SPUTTERED PARTICLES(',I1,') = ',I7,5X,
& 'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV')
1568 continue
do 1586 J=NJ(1)+1,JT(3)
write(21,1576) J-NJ(1),ISPIP(J),RIP(J),RIPJ(J),ESPIP(J),
& REIP(J),REIPJ(J),ESPMIP(J),J-NJ(1),ISPIS(J),RIS(J),
& RISJ(J),ESPIS(J),REIS(J),REISJ(J),ESPMIS(J),J-NJ(1),
& ISPOP(J),ROP(J),ROPJ(J),ESPOP(J),REOP(J),REOPJ(J),
& ESPMOP(J),J-NJ(1),ISPOS(J),ROS(J),ROSJ(J),ESPOS(J),
& REOS(J),REOSJ(J),ESPMOS(J)
1586 continue
write(21,1577)
do J=NJ(1)+1,JT(3)
write(21,1582) J-NJ(1),SPY(J),J,SPE(J),J,REY(J),J,EMSP(J)
enddo
1593 continue
do 7262 J=NJ(1)+1,JT(3)
if(IBSP(J).le.1) go to 7262
YSP=IBSP(J)
call MOMENTS(FIES0,SEES,THES,FOES,FIES,SIES,SIGMES,
& DFIES0,DSEES,DTHES,EBSP(J),SPE2S(J),SPE3S(J),
& SPE4S(J),SPE5S(J),SPE6S(J),YSP)
write(21,7117)
write(21,7261) J,FIES0,SEES,THES,FOES,SIGMES,DFIES0,
& DSEES,DTHES
FIESB(J)=FIES0
SEESB(J)=SEES
THESB(J)=THES
FOESB(J)=FOES
SGMESB(J)=SIGMES
DFIESB(J)=DFIES0
DSEESB(J)=DSEES
DTHESB(J)=DTHES
7262 continue
if(L.eq.2) go to 1532
if(ISPAL(3).eq.0) go to 1532
write(21,1534)
1534 format(//1X,'3.LAYER')
do 1536 J=JT(3)+1,LJ
write(21,1538) J-JT(3),IBSP(J),J-JT(3),EBSP(J)
1538 format(10X,'SPUTTERED PARTICLES(',I1,') = ',I7,6X,
& 'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV')
1536 continue
do J=JT(3)+1,LJ
write(21,1576) J-JT(3),ISPIP(J),RIP(J),RIPJ(J),ESPIP(J),
& REIP(J),REIPJ(J),ESPMIP(J),J-JT(3),ISPIS(J),RIS(J),
& RISJ(J),ESPIS(J),REIS(J),REISJ(J),ESPMIS(J),J-JT(3),
& ISPOP(J),ROP(J),ROPJ(J),ESPOP(J),REOP(J),REOPJ(J),
& ESPMOP(J),J-JT(3),ISPOS(J),ROS(J),ROSJ(J),ESPOS(J),
& REOS(J),REOSJ(J),ESPMOS(J)
enddo
write(21,1577)
do J=JT(3)+1,LJ
write(21,1582) J-JT(3),SPY(J),J-JT(3),SPE(J),
& J-JT(3),REY(J),J-JT(3),EMSP(J)
enddo
1532 continue
return
7117 format(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',
& 4X,' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',
& 3X,' ERROR 2.M ',3X,' ERROR 3.M ')
end
C WRITE TRANSMISSION SUMMARY
C=======================================================================
subroutine writeTransmissionSummary(nTrans,nPrimary,transEnergy,
& initialEnergy,elasticLossSum,elasticLossSquare,
& inElasticLossSum,inElasticLossSquare,energyMoment2,
& energyMoment3,energyMoment4,energyMoment5,energyMoment6,
& pathMoment1,pathMoment2,pathMoment3,pathMoment4,
& pathMoment5,pathMoment6,angleBins,angleX,
& angleDistribution,meanTransEnergy,transEnergySigma)
implicit none
integer nTrans
integer angleBins(20)
integer i
real*8 nPrimary,transEnergy,initialEnergy
real*8 elasticLossSum,elasticLossSquare
real*8 inElasticLossSum,inElasticLossSquare
real*8 energyMoment2,energyMoment3,energyMoment4
real*8 energyMoment5,energyMoment6
real*8 pathMoment1,pathMoment2,pathMoment3
real*8 pathMoment4,pathMoment5,pathMoment6
real*8 angleX(20),angleDistribution(20)
real*8 meanTransEnergy,transEnergySigma
real*8 transCount,particleTransmission,energyTransmission
real*8 relativeMeanEnergy,meanEnergy
real*8 meanElasticLoss,varElasticLoss,sigmaElasticLoss
real*8 errorElasticLoss
real*8 meanInElasticLoss,varInElasticLoss,sigmaInElasticLoss
real*8 errorInElasticLoss
real*8 energyVariance,energySkew,energyKurtosis
real*8 energyFifth,energySixth
real*8 energyError1,energyError2,energyError3
real*8 pathMean,pathVariance,pathSkew,pathKurtosis
real*8 pathFifth,pathSixth,pathSigma
real*8 pathError1,pathError2,pathError3
logical equal
if(nTrans.eq.0) return
write(21,1529)
1529 format(///5X,' TRANSMISSION OF PROJECTILES'/)
transCount=dble(nTrans)
particleTransmission=transCount/nPrimary
energyTransmission=transEnergy/(nPrimary*initialEnergy)
relativeMeanEnergy=energyTransmission/particleTransmission
meanEnergy=relativeMeanEnergy*initialEnergy
meanTransEnergy=meanEnergy
meanElasticLoss=0.0d0
varElasticLoss=0.0d0
sigmaElasticLoss=0.0d0
errorElasticLoss=0.0d0
meanInElasticLoss=0.0d0
varInElasticLoss=0.0d0
sigmaInElasticLoss=0.0d0
errorInElasticLoss=0.0d0
if (equal(transCount,1.0d0)) go to 1520
meanElasticLoss=elasticLossSum/transCount
varElasticLoss=elasticLossSquare/transCount
& -meanElasticLoss*meanElasticLoss
sigmaElasticLoss=dsqrt(varElasticLoss)
errorElasticLoss=sigmaElasticLoss/transCount
meanInElasticLoss=inElasticLossSum/transCount
varInElasticLoss=inElasticLossSquare/transCount
& -meanInElasticLoss*meanInElasticLoss
sigmaInElasticLoss=dsqrt(varInElasticLoss)
errorInElasticLoss=sigmaInElasticLoss/transCount
1520 write(21,1522) particleTransmission,energyTransmission,
& relativeMeanEnergy,meanEnergy
1522 format(//5X,'PART.TRANSM.COEF.=',1PE11.4,
& ' ENERGY TRANSM.COEF.=',1E11.4,
& ' REL.MEAN ENERGY =',1E11.4,
& ' MEAN ENERGY =',1E11.4)
call MOMENTS(meanTransEnergy,energyVariance,energySkew,
& energyKurtosis,energyFifth,energySixth,transEnergySigma,
& energyError1,energyError2,energyError3,transEnergy,
& energyMoment2,energyMoment3,energyMoment4,energyMoment5,
& energyMoment6,transCount)
call MOMENTS(pathMean,pathVariance,pathSkew,pathKurtosis,
& pathFifth,pathSixth,pathSigma,pathError1,pathError2,
& pathError3,pathMoment1,pathMoment2,pathMoment3,pathMoment4,
& pathMoment5,pathMoment6,transCount)
write(21,7117)
write(21,7241) meanTransEnergy,energyVariance,energySkew,
& energyKurtosis,transEnergySigma,energyError1,
& energyError2,energyError3
7241 format(1X,' ENERGY',5X,1P1E12.4,7E14.4)
write(21,7242) pathMean,pathVariance,pathSkew,pathKurtosis,
& pathSigma,pathError1,pathError2,pathError3
7242 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4)
write(21,7237) meanElasticLoss,varElasticLoss,sigmaElasticLoss,
& errorElasticLoss
write(21,7238) meanInElasticLoss,varInElasticLoss,
& sigmaInElasticLoss,errorInElasticLoss
write(21,1526)
1526 format(//5X,'POLAR ANGULAR DISTRIBUTION OF TRANSMITTED ',
& 'PARTICLES'//)
do i=1,20
angleDistribution(i)=dble(angleBins(i))*20.d0/dble(nTrans)
enddo
write(21,1530) (angleX(i),i=1,20),(angleBins(i),i=1,20),
& (angleDistribution(i),i=1,20)
1530 format(5X,20F6.2//,5X,20I6/5X,20F6.3)
return
7117 format(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',
& 4X,' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',
& 3X,' ERROR 2.M ',3X,' ERROR 3.M ')
7237 format(1X,'ELASTIC LOSS',5X,1P1E12.4,1E14.4,28X,
& 2E14.4)
7238 format(1X,' INEL. LOSS',5X,1P1E12.4,1E14.4,28X,
& 2E14.4)
end
C=======================================================================
C WRITE BACKSCATTERING SUMMARY
C=======================================================================
subroutine writeBackscatterSummary(nBackscatter,nBackLayer,
& nPrimary,backEnergySum,energyUnit,elasticLossSum,
& elasticLossSquare,inElasticLossSum,inElasticLossSquare,
& energyMoment2,energyMoment3,energyMoment4,energyMoment5,
& energyMoment6,logEnergyMoment1,logEnergyMoment2,
& logEnergyMoment3,logEnergyMoment4,logEnergyMoment5,
& logEnergyMoment6,pathMoment1,pathMoment2,pathMoment3,
& pathMoment4,pathMoment5,pathMoment6,angleBins,angleX,
& angleDistribution,meanBackEnergy,backEnergySigma)
implicit none
integer nBackscatter,nBackLayer
integer angleBins(20)
integer i
real*8 nPrimary,backEnergySum,energyUnit
real*8 elasticLossSum,elasticLossSquare
real*8 inElasticLossSum,inElasticLossSquare
real*8 energyMoment2,energyMoment3,energyMoment4
real*8 energyMoment5,energyMoment6
real*8 logEnergyMoment1,logEnergyMoment2,logEnergyMoment3
real*8 logEnergyMoment4,logEnergyMoment5,logEnergyMoment6
real*8 pathMoment1,pathMoment2,pathMoment3
real*8 pathMoment4,pathMoment5,pathMoment6
real*8 angleX(20),angleDistribution(20)
real*8 meanBackEnergy,backEnergySigma
real*8 bi,bil,particleReflection,energyReflection
real*8 relativeMeanEnergy,meanEnergy
real*8 meanElasticLoss,varElasticLoss,sigmaElasticLoss
real*8 errorElasticLoss
real*8 meanInElasticLoss,varInElasticLoss,sigmaInElasticLoss
real*8 errorInElasticLoss
real*8 energyVariance,energySkew,energyKurtosis
real*8 energyFifth,energySixth
real*8 energyError1,energyError2,energyError3
real*8 pathMean,pathVariance,pathSkew,pathKurtosis
real*8 pathFifth,pathSixth,pathSigma
real*8 pathError1,pathError2,pathError3
real*8 energyMoment1Out,energyMoment2Out,energyMoment3Out
real*8 energyMoment4Out,energyMoment5Out,energyMoment6Out
real*8 logMoment1Out,logMoment2Out,logMoment3Out
real*8 logMoment4Out,logMoment5Out,logMoment6Out
real*8 pathMoment1Out,pathMoment2Out,pathMoment3Out
real*8 pathMoment4Out,pathMoment5Out,pathMoment6Out
logical equal
if(nBackscatter.eq.0) return
write(21,1527)
1527 format(1H1,//5X,'BACKSCATTERING OF PROJECTILES'/)
bi=dble(nBackscatter)
bil=dble(nBackLayer)
particleReflection=bi/nPrimary
energyReflection=backEnergySum/(nPrimary*energyUnit)
relativeMeanEnergy=energyReflection/particleReflection
meanEnergy=backEnergySum/bi
meanBackEnergy=meanEnergy
if (equal(bi,1.0d0)) go to 1506
meanElasticLoss=elasticLossSum/bi
varElasticLoss=elasticLossSquare/bi
& -meanElasticLoss*meanElasticLoss
sigmaElasticLoss=dsqrt(varElasticLoss)
errorElasticLoss=sigmaElasticLoss/bi
meanInElasticLoss=inElasticLossSum/bi
varInElasticLoss=inElasticLossSquare/bi
& -meanInElasticLoss*meanInElasticLoss
sigmaInElasticLoss=dsqrt(varInElasticLoss)
errorInElasticLoss=sigmaInElasticLoss/bi
1506 write(21,1508) particleReflection,energyReflection,
& relativeMeanEnergy,meanEnergy
1508 format(/5X,'PART.REFL.COEF.=',1PE11.4,
& ' ENERGY REFL.COEF.=',1E11.4,
& ' REL.MEAN ENERGY =',1E11.4,
& ' MEAN ENERGY =',1E11.4)
call MOMENT(energyMoment1Out,energyMoment2Out,energyMoment3Out,
& energyMoment4Out,energyMoment5Out,energyMoment6Out,
& backEnergySum,energyMoment2,energyMoment3,energyMoment4,
& energyMoment5,energyMoment6,bi)
call MOMENT(logMoment1Out,logMoment2Out,logMoment3Out,
& logMoment4Out,logMoment5Out,logMoment6Out,
& logEnergyMoment1,logEnergyMoment2,logEnergyMoment3,
& logEnergyMoment4,logEnergyMoment5,logEnergyMoment6,bil)
call MOMENTS(meanBackEnergy,energyVariance,energySkew,
& energyKurtosis,energyFifth,energySixth,backEnergySigma,
& energyError1,energyError2,energyError3,backEnergySum,
& energyMoment2,energyMoment3,energyMoment4,energyMoment5,
& energyMoment6,bi)
call MOMENT(pathMoment1Out,pathMoment2Out,pathMoment3Out,
& pathMoment4Out,pathMoment5Out,pathMoment6Out,pathMoment1,
& pathMoment2,pathMoment3,pathMoment4,pathMoment5,
& pathMoment6,bi)
call MOMENTS(pathMean,pathVariance,pathSkew,pathKurtosis,
& pathFifth,pathSixth,pathSigma,pathError1,pathError2,
& pathError3,pathMoment1,pathMoment2,pathMoment3,pathMoment4,
& pathMoment5,pathMoment6,bi)
write(21,7117)
write(21,7241) meanBackEnergy,energyVariance,energySkew,
& energyKurtosis,backEnergySigma,energyError1,energyError2,
& energyError3
7241 format(1X,' ENERGY',5X,1P1E12.4,7E14.4)
write(21,7242) pathMean,pathVariance,pathSkew,pathKurtosis,
& pathSigma,pathError1,pathError2,pathError3
7242 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4)
write(21,7237) meanElasticLoss,varElasticLoss,sigmaElasticLoss,
& errorElasticLoss
write(21,7238) meanInElasticLoss,varInElasticLoss,
& sigmaInElasticLoss,errorInElasticLoss
write(21,7118)
write(21,1541) energyMoment1Out,energyMoment2Out,
& energyMoment3Out,energyMoment4Out,energyMoment5Out,
& energyMoment6Out
1541 format(1X,' ENERGY',5X,1P1E12.4,5E14.4)
write(21,1543) logMoment1Out,logMoment2Out,logMoment3Out,
& logMoment4Out,logMoment5Out,logMoment6Out
1543 format(1X,' LOGENERGY',5X,1P1E12.4,5E14.4)
write(21,1545) pathMoment1Out,pathMoment2Out,pathMoment3Out,
& pathMoment4Out,pathMoment5Out,pathMoment6Out
1545 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4)
do i=1,20
angleX(i)=0.05d0*dble(i)
enddo
write(21,1514)
1514 format(//5X,'POLAR ANGULAR DISTRIBUTION OF BACKSCATTERED ',
& 'PROJECTILES'//)
do i=1,20
angleDistribution(i)=dble(angleBins(i))*20.d0
& /dble(nBackscatter)
enddo
write(21,1518)(angleX(i),i=1,20),(angleBins(i),i=1,20),
& (angleDistribution(i),i=1,20)
1518 format(5X,20F6.2//,5X,20I6/5X,20F6.3)
return
7117 format(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',
& 4X,' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',
& 3X,' ERROR 2.M ',3X,' ERROR 3.M ')
7118 format(/20X,' 1.MOMENT ',4X,' 2.MOMENT ',4X,' 3.MOMENT ',
& 4X,' 4.MOMENT ',4X,' 5.MOMENT ',4X,' 6.MOMENT ')
7237 format(1X,'ELASTIC LOSS',5X,1P1E12.4,1E14.4,28X,
& 2E14.4)
7238 format(1X,' INEL. LOSS',5X,1P1E12.4,1E14.4,28X,
& 2E14.4)
end
C=======================================================================
C WRITE IMPLANTATION PROFILE AND RECOIL DEPTH TABLES
C=======================================================================
subroutine writeImplantationProfile(depth_interval_flag,YH,
& MAXD,MAXD1,LJ,NJ,L,XX,CW,IRL,IRP,RIRP,IPL,ION,DENT,
& DMGN,ELGD,PHON,CASMOT,ICD,ICDT,ICDR,ICDTR,ICDIRI,
& ICDIRN,IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,
& TCASMO,ICDTT,TIONR,TDENTR,TELGDR,TDMGNR,TPHONR,ICDTTR,
& ELE,ELI,ELP,ELD,ELET,ELIT,ELPT,ELDT,ICDJT,ICDJTR,
& ICDITR,ILD,DLI,IONR,DENTR,DMGNR,ELGDR,PHONR,MAXNL,
& MAXNL5)
implicit none
integer MAXD,MAXD1,LJ,L,IRL,MAXNL,MAXNL5
integer depth_interval_flag
integer NJ(MAXNL),IRP(0:MAXD1),IPL(MAXD)
integer ICD(MAXD,MAXNL5),ICDT(MAXD),ICDR(MAXD,MAXNL5)
integer ICDTR(MAXD),ICDIRI(MAXD,MAXNL5,MAXNL5)
integer ICDIRN(MAXD,MAXNL5),ICDJT(MAXNL5),ICDJTR(MAXNL5)
integer ICDITR(MAXNL5),ILD(MAXNL)
integer IIRP,IIPL,ICDTT,ICDTTR
integer I,J,K,im1,im2,im3
real*8 YH,CW,XX(MAXNL)
real*8 RIRP(0:MAXD1),ION(MAXD),DENT(MAXD),DMGN(MAXD)
real*8 ELGD(MAXD),PHON(MAXD),CASMOT(MAXD)
real*8 TION,TDENT,TDMGN,TELGD,TPHON,TCASMO,TRIRP
real*8 IONR(MAXD),DENTR(MAXD),DMGNR(MAXD),ELGDR(MAXD)
real*8 PHONR(MAXD),TIONR,TDENTR,TDMGNR,TELGDR,TPHONR
real*8 ELE(MAXD,MAXNL5),ELI(MAXD,MAXNL5)
real*8 ELP(MAXD,MAXNL5),ELD(MAXD,MAXNL5)
real*8 ELET(MAXNL5),ELIT(MAXNL5),ELPT(MAXNL5),ELDT(MAXNL5)
real*8 DLI(MAXNL)
real*8 D1,D2,Dmid
logical EQUAL
590 WRITE(21,600)
600 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES,2X,10HNORM.DEPTH,1X,
& 10HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X
& ,10HELLOSS(EV), 2X,10HDAMAGE(EV),2X,10HPHONON(EV),2X
& ,10HCASCAD(EV),5X,3HDPA/)
IF(depth_interval_flag.EQ.0) THEN
C WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ',
C & 'AGREE WITH LAYER THICKNESS'
WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT ',
& 'AGREE WITH LAYER THICKNESS'
ENDIF
#if defined (OS_WIN)
WRITE(22,6002)
6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES)
#else
write(22,'(a)') ' DEPTH PARTICLES'
#endif
IF(YH.LT.1.D0) GO TO 603
DO I=0,MAXD1
C This normalization of the implantation profile is wrong if the
C simulation does not account for all implanted projectiles,
C e.g. when the number of points in the profile is not enough to
C reach the maximum implantation depth.
RIRP(I) = DBLE(IRP(I))/YH
C To resolve this issue we need to check whether we reached maximum
C implantation or not. If not rerun calculation with larger step
C size!
ENDDO
603 D1=0.
D2=CW
WRITE(21,601) D1,IRP(0),RIRP(0)
601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4)
DO J=1,LJ
DO I=1,MAXD
ICDT(I)=ICDT(I)+ICD(I,J)
ICDTR(I)=ICDTR(I)+ICDR(I,J)
ENDDO
ENDDO
DO K=1,NJ(1)
DO J=1,NJ(1)
DO I=1,MAXD
ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J)
ENDDO
ENDDO
ENDDO
DO I=0,MAXD1
IIRP=IIRP+IRP(I)
TRIRP=TRIRP+RIRP(I)
ENDDO
DO I=1,MAXD
IIPL=IIPL+IPL(I)
TION=TION+ION(I)
TDENT=TDENT+DENT(I)
TDMGN=TDMGN+DMGN(I)
TELGD=TELGD+ELGD(I)
TPHON=TPHON+PHON(I)
TCASMO=TCASMO+CASMOT(I)
ICDTT=ICDTT+ICDT(I)
TIONR=TIONR+IONR(I)
TDENTR=TDENTR+DENTR(I)
TELGDR=TELGDR+ELGDR(I)
TDMGNR=TDMGNR+DMGNR(I)
TPHONR=TPHONR+PHONR(I)
ICDTTR=ICDTTR+ICDTR(I)
ENDDO
do im1=MAXD,1,-1
if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20
enddo
im1=1
20 im1=min0(im1+2,MAXD)
DO I=1,im1
WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I),
& DMGN(I),ELGD(I),PHON(I),CASMOT(I),ICDT(I)
Dmid=(D2-D1)/2+D1
WRITE(22,701) Dmid,IRP(I)
700 FORMAT(1X,F6.0,1H-,F6.0,I10,E12.4,I10,1P1E14.4,5E12.4,I8)
701 FORMAT(1X, F16.4, 2X, I10)
D1=D2
D2=D2+CW
ENDDO
WRITE(21,604) D2-CW,IRP(MAXD1),RIRP(MAXD1)
604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4)
WRITE(21,710) IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO
& ,ICDTT
710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8)
DO J=1,NJ(1)
DO I=1,MAXD
ELET(J)=ELET(J)+ELE(I,J)
ELIT(J)=ELIT(J)+ELI(I,J)
ELPT(J)=ELPT(J)+ELP(I,J)
ELDT(J)=ELDT(J)+ELD(I,J)
ICDJT(J)=ICDJT(J)+ICD(I,J)
ICDJTR(J)=ICDJTR(J)+ICDR(I,J)
ICDITR(J)=ICDITR(J)+ICDIRN(I,J)
ENDDO
ENDDO
WRITE(21,1521)
1521 FORMAT(1H1,4X,'DEPTH(A)'
& ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)'
& ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)'
& ,2X,'DPA(1)',2X,'DPA(2)'/)
D1=0.
D2=CW
do im2=MAXD,1,-1
if(.NOT.EQUAL(eli(im2,1),0.D0).or.(.NOT.EQUAL(eli(im2,2),0.D0))
& ) goto 30
enddo
im2=1
30 im2=MIN0(im2+2,MAXD)
DO 1525 I=1,im2
WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1),
& ELI(I,2),ELE(I,2),ELD(I,2),ELP(I,2),ICD(I,1),ICD(I,2)
1523 FORMAT(1X,F6.0,1H-,F6.0,1P8E12.4,2I8)
D1=D2
D2=D2+CW
1525 CONTINUE
WRITE(21,1533) ELIT(1),ELET(1),ELDT(1),ELPT(1),ELIT(2),ELET(2)
& ,ELDT(2),ELPT(2),ICDJT(1),ICDJT(2)
1533 FORMAT(/14X,1P8E12.4,2I8///)
DO I=1,L-1
ILD(I)=IDINT(XX(I)/CW+0.01D0)
IF(ILD(I).GT.MAXD) ILD(I)=MAXD
DO J=1,ILD(I)
DLI(I)=DLI(I)+DMGN(J)
ENDDO
ENDDO
DLI(L)=TDMGN
DO 1493 I=L,2,-1
DLI(I)=DLI(I)-DLI(I-1)
1493 CONTINUE
DO 1494 I=1,L
WRITE(21,1495) I,DLI(I)
1495 FORMAT(/5X,'DAMAGE IN LAYER ',I3,' : ',1P1E12.4)
1494 CONTINUE
1455 CONTINUE
if(irl.eq.0) goto 1497
WRITE(21,1496)
1496 FORMAT(1H1,/,5X,'RECOILS')
WRITE(21,1597)
1597 FORMAT(///,5X,8HDEPTH(A), 5X,10HINLOSS(EV),3X,10HTELOSS(EV),3X
& ,10HELLOSS(EV), 3X,10HDAMAGE(EV),3X,10HPHONON(EV),5X,3HDPA,
& 2X,6HDPA(1),2X,6HDPA(2), 1X,5H(1-1),1X,5H(1-2),1X,5H(2-1),1X
& ,5H(2-2)/)
D1=0.D0
D2=CW
do im3=MAXD,1,-1
if (.not.equal(ionr(im3),0.D0)) go to 31
enddo
im3=1
31 im3=MIN0(im3+2,MAXD)
DO I=1,im3
WRITE(21,1595) D1,D2,IONR(I),DENTR(I),DMGNR(I),ELGDR(I),
& PHONR(I),ICDTR(I),ICDIRN(I,1),ICDIRN(I,2)
& ,ICDIRI(I,1,1),ICDIRI(I,1,2),ICDIRI(I,2,1),ICDIRI(I,2,2)
1595 FORMAT(1X,F6.0,1H-,F6.0,1P1E14.4,4E13.4,3I8,4I6)
D1=D2
D2=D2+CW
ENDDO
WRITE(21,1596) TIONR,TDENTR,TDMGNR,TELGDR,TPHONR
& ,ICDTTR,ICDITR(1),ICDITR(2)
1596 FORMAT(/14X,1P1E14.4,4E13.4,3I8)
1497 continue
return
end
C=======================================================================
C WRITE INPUT AND BEAM SUMMARY
C=======================================================================
subroutine writeInputSummary(inputName,outputName,rangeName,
& zProjectile,massProjectile,energy0,energySigma,alpha,
& alphaSigma,finalEnergy,surfaceEnergy,sheath,recoilCutoff,
& nProjectiles,ri,ri2,ri3,xStart,rd,depthBin,ca,kk0,kk0r,
& stoppingPrimary,stoppingRecoil,potentialPrimary,
& potentialRecoil,recoilMode,nLayers,layerThickness,density,
& correction,nElements,zTarget,massTarget,concentration,
& surfaceBinding,displacementEnergy,bulkBinding,ch1,ch2,ch3,
& ch4,ch5,ti,zArg,velc,emx,maxNl,maxEl,secondsStartTotal)
implicit none
integer nProjectiles,kk0,kk0r,stoppingPrimary,stoppingRecoil
integer potentialPrimary,potentialRecoil,recoilMode,nLayers
integer maxNl,maxEl,secondsStartTotal,secondsStopTotal
integer nElements(maxNl)
integer potKrc,potMoliere,potZbl
integer stopLs,stopOr,stopMixed,stopIcru49,stopZiegler
integer i,j
real*8 zProjectile,massProjectile,energy0,energySigma,alpha
real*8 alphaSigma,finalEnergy,surfaceEnergy
& ,sheath,recoilCutoff
real*8 ri,ri2,ri3,xStart,rd,depthBin,ca
real*8 layerThickness(maxNl),density(maxNl),correction(maxNl)
real*8 zTarget(maxNl,maxEl),massTarget(maxNl,maxEl)
real*8 concentration(maxNl,maxEl)
real*8 surfaceBinding(maxNl,maxEl)
real*8 displacementEnergy(maxNl,maxEl)
real*8 bulkBinding(maxNl,maxEl)
real*8 ch1(maxNl,maxEl),ch2(maxNl,maxEl)
real*8 ch3(maxNl,maxEl),ch4(maxNl,maxEl),ch5(maxNl,maxEl)
real*8 ti,zArg,velc,emx
character*(*) inputName,outputName,rangeName
character*18 potentialName,potentialRecoilName
character*18 stoppingName,stoppingRecoilName
character monthStop*4,dayStop*2,yearStop*4,hourStop*2
character minStop*2,secStop*2
logical EQUAL
parameter (potKrc=1,potMoliere=2,potZbl=3)
parameter (stopLs=1,stopOr=2,stopMixed=3)
parameter (stopIcru49=4,stopZiegler=5)
C How many seconds are needed for the simulation ??
CALL TimeStamp(dayStop,monthStop,yearStop,hourStop,minStop,
& secStop,secondsStopTotal)
WRITE(21,*)
WRITE(21,10051)dayStop,monthStop,yearStop,hourStop,minStop,
& secStop
WRITE(*,10051)dayStop,monthStop,yearStop,hourStop,minStop,
& secStop
10051 FORMAT(1x,' End: ',A2,'.',A4,1x,A4,1x,A2,':',A2,':',A2)
WRITE(21,*)
WRITE(21,10052)nProjectiles,(secondsStopTotal-secondsStartTotal)
10052 FORMAT(1x,' Simulation needed for ',I7,' projectiles ',I7,
& ' seconds')
WRITE(21,1402)inputName
1402 FORMAT(//30X,'* INPUT DATA *',5X,A12)
WRITE(21,1404) zProjectile,massProjectile,energy0,
& energySigma,alpha,alphaSigma,finalEnergy,surfaceEnergy,
& sheath,recoilCutoff
1404 FORMAT(//,7X,2HZ1,8X,2HM1,10X,2HE0,6X,4HEsig,7X,5HALPHA,7X
& ,8HALPHASIG,7X,2HEF,7X ,3HESB,6X,6HSHEATH,5X,3HERC/2F10.2
& ,1F13.2,7F10.2)
WRITE(21,1406) nProjectiles,ri,ri2,ri3,xStart,rd,depthBin,ca,
& kk0,kk0r,stoppingPrimary,stoppingRecoil,potentialPrimary,
& potentialRecoil,recoilMode
1406 FORMAT(/7X,2HNH,8X,2HRI,5X,3HRI2,5X,3HRI3,11X,2HX0,8X,2HRD,8X,2HCW
& ,8X,2HCA ,7X,3HKK0,3X,4HKK0R,3X,5HKDEE1,2X,5HKDEE2,2X,4HIPOT
& ,3X,5HIPOTR ,3X,3HIRL/I10,3F10.2,1F13.2,3F10.2,1X,7I7)
WRITE(21,1408)
1408 FORMAT(//13X,2HDX,6X,3HRHO,4X,2HCK,2X
& ,5HZ(,1),1X,5HZ(,2),1X,5HZ(,3),1X,5HZ(,4),1X,5HZ(,5),2X
& ,5HM(,1),2X,5HM(,2),2X,5HM(,3),2X,5HM(,4),2X,5HM(,5),1X
& ,5HC(,1),1X,5HC(,2),1X,5HC(,3),1X,5HC(,4),1X,5HC(,5))
DO I=1,nLayers
WRITE(21,1412) I,layerThickness(I),density(I),correction(I),
& (zTarget(I,J),J=1,5),(massTarget(I,J),J=1,5),
& (concentration(I,J),J=1,5)
IF(nElements(I).GT.5) THEN
WRITE(21,5113)
DO J=6,nElements(I)
WRITE(21,5111) J,zTarget(I,J),massTarget(I,J),
& concentration(I,J)
ENDDO
ENDIF
1412 FORMAT(/1X,I3,6H.LAYER,1X,1F8.2,2F7.2,5F6.0,5F7.2,5F6.3)
5113 FORMAT(6X,'ADDITIONAL ELEMENTS: J Z M C')
5111 FORMAT(6X,23X,I3,1X,F6.0,1X,F7.2,1X,F8.5)
5114 FORMAT(6X,'ADDITIONAL ELEMENTS: J SBE ED BE')
5115 FORMAT(6X,23X,I3,1X,F6.2,1X,F7.2,1X,F8.2)
ENDDO
WRITE(21,1414)
1414 FORMAT(//27X,'***',2X,'SBE(LAYER,ELEMENT)',2X,'***',5X
& ,'***',5X,'ED(LAYER,ELEMENT)',5X,'***',5X
& ,'***',3X,'BE(LAYER,ELEMENT)',2X,'***')
DO 1416 I=1,nLayers
WRITE(21,1418) I,(surfaceBinding(I,J),J=1,5),
& (displacementEnergy(I,J),J=1,5),(bulkBinding(I,J),J=1,5)
IF(nElements(I).GT.5) THEN
WRITE(21,5114)
DO J=6,nElements(I)
WRITE(21,5115) J,surfaceBinding(I,J),
& displacementEnergy(I,J),bulkBinding(I,J)
ENDDO
ENDIF
1418 FORMAT(/1X,I3,6H.LAYER,17X,5F6.2,3X,5F7.2,3X,5F6.2)
1416 CONTINUE
IF(stoppingPrimary.LT.stopIcru49) GO TO 1421
WRITE(21,1419)
1419 FORMAT(//30X,'CH1',10X,'CH2',10X,'CH3',10X,'CH4',10X,'CH5')
DO 1417 I=1,nLayers
DO J=1,nElements(I)
IF(J.EQ.1) THEN
WRITE(21,1415) I,ch1(I,J),ch2(I,J),ch3(I,J),
& ch4(I,J),ch5(I,J)
ELSE
WRITE(21,1423) ch1(I,J),ch2(I,J),ch3(I,J),
& ch4(I,J),ch5(I,J)
ENDIF
ENDDO
1417 CONTINUE
1415 FORMAT(/1X,I3,6H.LAYER,17X,5F13.6)
1423 FORMAT(/25X,5F13.6)
1421 CONTINUE
IF(potentialPrimary.EQ.potKrc) potentialName='KR-C POTENTIAL'
IF(potentialPrimary.EQ.potMoliere) potentialName='mod. MOLIERE '
IF(potentialPrimary.EQ.potZbl) potentialName='ZBL POTENTIAL'
IF(potentialRecoil.EQ.potKrc) potentialRecoilName='KR-C POTENTIAL'
IF(potentialRecoil.EQ.potMoliere) potentialRecoilName=
& 'MOLIERE POTENTIAL'
IF(potentialRecoil.EQ.potZbl) potentialRecoilName='ZBL POTENTIAL'
WRITE(21,1411) potentialName,potentialRecoilName
1411 FORMAT(//7X,'INTERACTION POTENTIAL : PROJECTILE-TARGET : ',A18
& ,' TARGET-TARGET : ',A18)
IF(stoppingPrimary.EQ.stopLs) stoppingName='LINDHARD-SCHARFF'
IF(stoppingPrimary.EQ.stopOr) stoppingName='OEN-ROBINSON'
IF(stoppingPrimary.EQ.stopMixed) stoppingName='50% LS 50% OR'
IF(stoppingPrimary.EQ.stopIcru49) stoppingName='AZ nach ICRU49'
IF(stoppingPrimary.EQ.stopZiegler) stoppingName='ZIEGLER'
IF(stoppingRecoil.EQ.stopLs) stoppingRecoilName='LINDHARD-SCHARFF'
IF(stoppingRecoil.EQ.stopOr) stoppingRecoilName='OEN-ROBINSON'
IF(stoppingRecoil.EQ.stopMixed) stoppingRecoilName='50% LS 50% OR'
WRITE(21,1413) stoppingName,stoppingRecoilName
1413 FORMAT(//7X,'INELASTIC LOSS MODEL : PROJECTILE-TARGET : ',A18
& ,' TARGET-TARGET : ',A18)
IF(energy0.GT.0.D0) GO TO 1420
IF(alpha.LT.0.D0) GO TO 1405
WRITE(21,1422) ti,zArg,velc,emx
1422 FORMAT(//6X,'MAXWELLIAN DISTRIBUTION',7X,2HTI,5X,4HZARG,5X ,4HVELC
& ,8X,3HEMX/29X,1F10.2,2F9.4,1E14.6)
GO TO 1427
1405 continue
WRITE(21,1407) ti,sheath,-alpha,emx
1407 FORMAT(//6X,'MAXWELLIAN DISTRIBUTION (ENERGY)',7X,'TI',5X
& ,'SHEATH',5X,'ALPHAM',8X,'EMX'/38X,3F10.2,2X,1E14.6)
GO TO 1427
1420 IF(alpha.EQ.-1.) WRITE(21,1424)
1424 FORMAT(//6X,'RANDOM DISTRIBUTION'/)
IF(alpha.EQ.-2.) WRITE(21,1426)
1426 FORMAT(//6X,'COSINE DISTRIBUTION'/)
1427 CONTINUE
IF(EQUAL(energySigma,0.D0)) THEN
WRITE(21,14271)
14271 FORMAT(//6X,'fixed PROJECTILE ENERGY'/)
ELSE
WRITE(21,14272)
14272 FORMAT(//6X,'PROJECTILE ENERGY has GAUSSIAN DISTRIBUTION '/)
ENDIF
IF(EQUAL(alphaSigma,0.D0)) THEN
WRITE(21,14273)
14273 FORMAT(//6X,'fixed PROJECTILE ANGLE'/)
ELSE
WRITE(21,14274)
14274 FORMAT(//6X,'PROJECTILE ANGLE has 1D GAUSSIAN DISTRIBUTION '/)
ENDIF
WRITE(21,1428)outputName
1428 FORMAT(1H1,//30X,'* OUTPUT DATA *',5X,A12)
#if defined (OS_WIN)
WRITE(22,14280)rangeName
14280 FORMAT(1H1,//30X,'* RANGE DATA *',5X,A12)
#endif
RETURN
END
C=======================================================================
C WRITE TARGET SETUP AND PRECOMPUTED TABLES
C=======================================================================
subroutine writeTargetSetupSummary(hlm,hlmt,su,sut,xc,rt,sfe,
& inelasticMode,nLayers,nSpecies,nProj,kib,kit,maxActive,nAll,
& nPrimary,nSecondary,kis,kist,layerEps0,layerZ,layerMass,
& layerArealDensity,layerLm,layerPdMax,layerASig,layerBinding,
& layerBoundary,nElements,a1,kor1,a,kor,maxNl,maxSpecies)
implicit none
integer maxNl,maxSpecies
integer*4 inelasticMode,nLayers,nSpecies,nProj,kib,kit
integer*4 maxActive,nAll,nPrimary,nSecondary,kis,kist
integer*4 nElements(maxNl)
integer i,j
real*8 hlm,hlmt,su,sut,xc,rt,sfe
real*8 layerEps0(maxNl),layerZ(maxNl),layerMass(maxNl)
real*8 layerArealDensity(maxNl),layerLm(maxNl)
real*8 layerPdMax(maxNl),layerASig(maxNl)
real*8 layerBinding(maxNl),layerBoundary(maxNl)
real*8 a1(maxSpecies),kor1(maxSpecies)
real*8 a(maxSpecies,maxSpecies),kor(maxSpecies,maxSpecies)
C
C This block reports the layer setup and precomputed pair tables used
C by the transport calculation. It is output-only and must not change
C the simulation state.
C
write(21,1430) hlm,hlmt,su,sut,xc,rt,sfe,inelasticMode,
& nLayers,nSpecies
1430 format(//17X,'HLM',7X,'HLMT',8X,'SU',7X,'SUT',8X,'XC',8X,
& 'RT',7X,'SFE',6X,'INEL',9X,'L',8X,'LJ'/10X,1F11.4,
& 1F10.3,1F10.4,1F10.3,1F10.4,1F10.3,1F10.2,3I10)
write(21,1432) nProj,kib,kit,maxActive,nAll,nPrimary,
& nSecondary,kis,kist
1432 format(//16X,'NPROJ',7X,'KIB',7X,'KIT',6X,'MAXA',6X,
& 'NALL',7X,'NPA',7X,'NSA',7X,'KIS',6X,'KIST'/11X,9I10)
write(21,470)
470 format(//12X,'EPS0(I)',7X,'Z2(I)',7X,'M2(I)',5X,
& 'ARHO(I)' ,7X,'LM(I)',5X,'PDMAX(I)',5X,'ASIG(I)',7X,
& 'SB(I)',7X,'XX(I)' ,8X,'NJ(I)')
do i=1,nLayers
write(21,473) i,layerEps0(i),layerZ(i),layerMass(i),
& layerArealDensity(i),layerLm(i),layerPdMax(i),
& layerASig(i),layerBinding(i),layerBoundary(i),
& nElements(i)
473 format(/1X,I3,6H.LAYER,1X,9E12.4,I10)
enddo
write(21,474)
474 format(//13X,
& 'A1(1)',3X,'A1(2)',3X,'A1(3)',3X,'A1(4)',3X,
& 'A1(5)',3X,'A1(6)',3X,'A1(7)',3X,'A1(8)',3X,
& 'A1(9)',2X,'A1(10)',2X,'A1(11)',2X,'A1(12)',2X,
& 'A1(13)',2X,'A1(14)',2X,'A1(15)',2X,'A1(16)',2X,
& 'A1(17)',2X,'A1(18)',2X,'A1(19)',2X,'A1(20)',2X,
& 'A1(21)',2X,'A1(22)',2X,'A1(23)',2X,'A1(24)',2X,
& 'A1(25)',2X,'A1(26)',2X,'A1(27)',2X,'A1(28)',2X,
& 'A1(29)',2X,'A1(30)',2X,'A1(31)',2X,'A1(32)',2X,
& 'A1(33)',2X,'A1(34)',2X,'A1(35)')
do 475 i=1,nSpecies
write(21,471) a1(i)
471 format(/1X,9X,35F8.5)
475 continue
write(21,484)
484 format(//11X,
& 'KOR1(1)',1X,'KOR1(2)',1X,'KOR1(3)',1X,'KOR1(4)',1X,
& 'KOR1(5)',1X,'KOR1(6)',1X,'KOR1(7)',1X,'KOR1(8)',1X,
& 'KOR1(9)',1X,'KOR1(A)',1X,'KOR1(B)',1X,'KOR1(C)',1X,
& 'KOR1(D)',1X,'KOR1(E)',1X,'KOR1(F)',1X,'KOR1(G)',1X,
& 'KOR1(H)',1X,'KOR1(I)',1X,'KOR1(J)',1X,'KOR1(K)',1X,
& 'KOR1(L)',1X,'KOR1(M)',1X,'KOR1(N)',1X,'KOR1(O)',1X,
& 'KOR1(P)',1X,'KOR1(Q)',1X,'KOR1(R)',1X,'KOR1(S)',1X,
& 'KOR1(T)',1X,'KOR1(U)',1X,'KOR1(V)',1X,'KOR1(W)',1X,
& 'KOR1(X)',1X,'KOR1(Y)',1X,'KOR1(Z)')
do 486 i=1,nSpecies
write(21,489) kor1(i)
489 format(/1X,9X,35F8.5)
486 continue
write(21,476)
476 format(//12X,
& 'A(I,1)',2X,'A(I,2)',2X,'A(I,3)',2X,'A(I,4)',2X,
& 'A(I,5)',2X,'A(I,6)',2X,'A(I,7)',2X,'A(I,8)',2X,
& 'A(I,9)',1X,'A(I,10)',1X,'A(I,11)',1X,'A(I,12)',1X,
& 'A(I,13)',1X,'A(I,14)',1X,'A(I,15)',1X,'A(I,16)',1X,
& 'A(I,17)',1X,'A(I,18)',1X,'A(I,19)',1X,'A(I,20)',1X,
& 'A(I,21)',1X,'A(I,22)',1X,'A(I,23)',1X,'A(I,24)',1X,
& 'A(I,25)',1X,'A(I,26)',1X,'A(I,27)',1X,'A(I,28)',1X,
& 'A(I,29)',1X,'A(I,30)',1X,'A(I,31)',1X,'A(I,32)',1X,
& 'A(I,33)',1X,'A(I,34)',1X,'A(I,35)')
do i=1,nSpecies
write(21,477) (a(i,j),j=1,nSpecies)
477 format(/1X,9X,35F8.5)
enddo
write(21,490)
490 format(//11X,
& 'KOR(,1)',1X,'KOR(,2)',1X,'KOR(,3)',1X,'KOR(,4)',1X,
& 'KOR(,5)',1X,'KOR(,6)',1X,'KOR(,7)',1X,'KOR(,8)',1X,
& 'KOR(,9)',1X,'KOR(,A)',1X,'KOR(,B)',1X,'KOR(,C)',1X,
& 'KOR(,D)',1X,'KOR(,E)',1X,'KOR(,F)',1X,'KOR(,G)',1X,
& 'KOR(,H)',1X,'KOR(,I)',1X,'KOR(,J)',1X,'KOR(,K)',1X,
& 'KOR(,L)',1X,'KOR(,M)',1X,'KOR(,N)',1X,'KOR(,O)',1X,
& 'KOR(,P)',1X,'KOR(,Q)',1X,'KOR(,R)',1X,'KOR(,S)',1X,
& 'KOR(,T)',1X,'KOR(,U)',1X,'KOR(,V)',1X,'KOR(,W)',1X,
& 'KOR(,X)',1X,'KOR(,Y)',1X,'KOR(,Z)')
do i=1,nSpecies
write(21,492) (kor(i,j),j=1,nSpecies)
492 format(/1X,9X,35F8.5)
enddo
return
end
C=======================================================================
C WRITE INTEGRAL PARTICLE SUMMARY
C=======================================================================
subroutine writeIntegralSummary(nImplanted,implantedEnergy,
& nBackscatter,backEnergy,nTransmitted,transEnergy,
& nBackSputter,backSputterEnergy,nTransSputter,
& transSputterEnergy,tryEnergy,negativeEnergy,reducedEnergy,
& reflectionCoeff,nLayers,layerBoundary,binWidth,
& implantedProfile,implantedLayerBins,layerBinEnd,
& numberInLayer,depthIntervalFlag,meanElasticCollisions,
& meanDisplacingCollisions,meanSurfaceCollisions,maxNl,
& maxD)
implicit none
integer nImplanted,nBackscatter,nTransmitted,nBackSputter
integer nTransSputter,tryEnergy,negativeEnergy,nLayers
integer depthIntervalFlag,maxNl,maxD,i,j
integer implantedProfile(0:maxD-1),implantedLayerBins(maxNl)
integer layerBinEnd(0:maxNl),numberInLayer(maxNl)
real*8 implantedEnergy,backEnergy,transEnergy,backSputterEnergy
real*8 transSputterEnergy,reducedEnergy,reflectionCoeff
real*8 layerBoundary(maxNl),binWidth
real*8 meanElasticCollisions,meanDisplacingCollisions
real*8 meanSurfaceCollisions
C
C Main integral counts and energies for implanted, backscattered,
C transmitted, and sputtered particles.
C
write(21,1500) nImplanted,implantedEnergy,nBackscatter,
& backEnergy,nTransmitted,transEnergy,nBackSputter,
& backSputterEnergy,nTransSputter,transSputterEnergy,
& tryEnergy,negativeEnergy,reducedEnergy,reflectionCoeff
1500 format(1H1,//11X,20HIMPLANTED PARTICLES=,I7,5X,7HENERGY=,
& E10.4,3H EV/7X,24HBACKSCATTERED PARTICLES=,I7,5X,
& 7HENERGY=,E10.4,3H EV/9X,22HTRANSMITTED PARTICLES=,
& I7,5X,7HENERGY=,E10.4,3H EV/7X,
& 24HBACKSPUTTERED PARTICLES=,I7,5X,7HENERGY=,E10.4,
& 3H EV/6X,'TRANSM. SPUTT. PARTICLES=',I7,5X,
& 7HENERGY=,E10.4,3H EV/15X,16HTRIED PARTICLES=,I7/
& 9X,22HPARTICLES with neg. E=,I7/6X,
& 25HTHOMAS FERMI RED. ENERGY=,2X,E10.4/6X,
& 25HSCALED PART. REFL. COEFF=,2x,E10.4)
if(nLayers.gt.1) then
do i=1,nLayers
layerBinEnd(i)=IDINT(layerBoundary(i)/binWidth+0.01)
enddo
do i=1,nLayers
do j=layerBinEnd(i-1)+1,layerBinEnd(i)
implantedLayerBins(i)=implantedLayerBins(i)
& + implantedProfile(j)
enddo
enddo
write(21,15001)
15001 format(/33x,'FROM BIN WIDTH',2x,'FROM LAYER THICKNESS')
if(depthIntervalFlag.eq.0) write(21,'(36x,5HWRONG)')
do i=1,nLayers
write(21,1501) i,implantedLayerBins(i),numberInLayer(i)
enddo
1501 format(/1x,'IMPLANTED PARTICLES (LAYER ',i1,')=',i16,
& 2x,i16)
endif
write(21,1511) meanElasticCollisions,meanDisplacingCollisions,
& meanSurfaceCollisions
1511 format(//2X,'PROJECTILES : ',
& 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X,
& 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F8.3/65X,
& 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F8.3)
return
end
C=======================================================================
C SET UP MATRIX TABLE AXES
C=======================================================================
subroutine setupMatrixTableAxes(energyBinWidth,angleBinWidth,
& polarBinWidth,nProjectiles,nPolarBinsMinusOne,
& nAngleBinsMinusOne,nAzimuthBins,maxDepthBin,jt,meab,magb,
& ema,meat,magt,emat,mepb,mept,mags,magst,magsa,meas,
& meast,measl,meastl,meags,meagb,meagt,energyLog,
& tempNh,maxNl,maxD2,maxSpeciesMinusFirstLayer)
implicit none
integer nProjectiles,nPolarBinsMinusOne,nAngleBinsMinusOne
integer nAzimuthBins,maxDepthBin,maxNl,maxD2
integer maxSpeciesMinusFirstLayer
integer j,k,azimuthBin
integer jt(maxNl)
integer meab(maxD2,22),magb(62,22),meagb(maxD2,36,22)
integer meat(maxD2,22),magt(62,22),meagt(maxD2,36,22)
integer mepb(maxD2,maxD2),mept(maxD2,maxD2)
integer meas(maxD2,22,maxSpeciesMinusFirstLayer)
integer meast(maxD2,22,maxSpeciesMinusFirstLayer)
integer mags(62,22,maxSpeciesMinusFirstLayer)
integer magst(62,22,maxSpeciesMinusFirstLayer)
integer magsa(62,32,maxSpeciesMinusFirstLayer)
integer measl(75,21,maxSpeciesMinusFirstLayer)
integer meastl(75,21,maxSpeciesMinusFirstLayer)
integer meags(maxD2,12,22,maxSpeciesMinusFirstLayer)
real*8 energyBinWidth,angleBinWidth,polarBinWidth
real*8 ema(62,22),emat(62,22),energyLog(75),tempNh
real*8 logBinWidth
C
C Fill the first row/column labels used by the matrix-style output
C tables. The simulation accumulates matrix entries during transport;
C this routine only sets the printable axes before those tables are
C written.
C
do j=2,nPolarBinsMinusOne
magb(j,1) = (j-1)*polarBinWidth
magt(j,1) = (j-1)*polarBinWidth
ema(j,1)=DBLE(j-1)*polarBinWidth
emat(j,1)=DBLE(j-1)*polarBinWidth
enddo
do j=2,21
meab(1,j) = j-1
meat(1,j) = j-1
magb(1,j) = j-1
magt(1,j) = j-1
ema(1,j) = j-1
emat(1,j) = j-1
enddo
do j=2,maxDepthBin
meab(j,1) = j-1
meat(j,1) = j-1
mepb(j,1) = j-1
mepb(1,j) = j-1
mept(j,1) = j-1
mept(1,j) = j-1
enddo
do k=1,jt(3)
do j=2,nPolarBinsMinusOne
mags(j,1,k) = (j-1)*polarBinWidth
magst(j,1,k) = (j-1)*polarBinWidth
magsa(j,1,k) = (j-1)*polarBinWidth
enddo
do j=2,nAngleBinsMinusOne
magsa(1,j,k) = (j-1)*angleBinWidth
enddo
do j=2,21
meas(1,j,k) = j-1
meast(1,j,k) = j-1
mags(1,j,k) = j-1
magst(1,j,k) = j-1
enddo
do j=2,maxDepthBin
meas(j,1,k) = j-1
meast(j,1,k) = j-1
enddo
do j=1,20
measl(1,j,k)=j
meastl(1,j,k)=j
enddo
do azimuthBin=1,nAzimuthBins,1
do j=2,21
meags(1,azimuthBin,j,k) = j-1
enddo
do j=2,maxDepthBin
meags(j,azimuthBin,1,k) = j-1
enddo
enddo
enddo
do azimuthBin=1,nAzimuthBins,1
do j=2,21
meagb(1,azimuthBin,j) = j-1
meagt(1,azimuthBin,j) = j-1
enddo
do j=2,maxDepthBin
meagb(j,azimuthBin,1) = j-1
meagt(j,azimuthBin,1) = j-1
enddo
enddo
do j=2,74
energyLog(j)=10.D0**(j/12.D0)*10.D0**(-7.D0/6.D0)
enddo
logBinWidth=(10.D0**(1.D0/12.D0)-1.D0)*10.D0**(-7.D0/6.D0)
tempNh=logBinWidth*DBLE(nProjectiles)
return
end
C=======================================================================
C WRITE GUI SEQUENCE SUMMARY FILE
C=======================================================================
subroutine writeSequenceSummary(seqNam,chem,nLayers,e0KeV,
& esigKeV,alpha,alphaSig,nProjectiles,nImplanted,
& nBackscatter,nTransmitted,tryEnergy,negativeEnergy,
& meanRange,rangeSigma,meanBackEnergy,backEnergySigma,
& meanTransEnergy,transEnergySigma,reducedEnergy,
& reflectionCoeff,numberInLayer,maxNl)
implicit none
integer maxNl,nLayers,nProjectiles,nImplanted,nBackscatter
integer nTransmitted,tryEnergy,negativeEnergy,k
integer numberInLayer(maxNl)
real*8 e0KeV,esigKeV,alpha,alphaSig
real*8 meanRange,rangeSigma,meanBackEnergy,backEnergySigma
real*8 meanTransEnergy,transEnergySigma,reducedEnergy
real*8 reflectionCoeff
character*64 seqNam
character*50 chem(maxNl)
C
C The *.seq file is the compact per-run summary consumed by the GUI.
C Keep this writer separate from the long human-readable output report.
C
open(unit=33,file=seqNam,status='replace')
write(33,7802) (chem(k),k=1,nLayers)
write(33,7801) e0KeV,esigKeV,alpha,alphaSig,nProjectiles,
& nImplanted,nBackscatter,nTransmitted,tryEnergy,
& negativeEnergy,meanRange,rangeSigma,meanBackEnergy,
& backEnergySigma,meanTransEnergy,transEnergySigma,
& reducedEnergy,reflectionCoeff,
& (numberInLayer(k),k=1,nLayers)
close(33)
return
7802 format(6x,'Energy',4x,'SigmaE',5x,'Alpha',2x,'SigAlpha',
& 4x,'ntot',5x,'imp',2x,'backsc',3x,'trans',3x,
& 'tried',4x,'negE',3x,'range',6x,'straggeling',2x,
& 'Eback',7x,'sigEback',4x,'Etrans',6x,'SigEtrans',3x,
& 'red._E',6x,'PRC',3x,999(A20,1x))
7801 format(F12.2,3(1x,F9.2),1x,6(I7,1x),6(E12.4),
& 2(E12.4),999(I7,1x))
end
C=======================================================================
C READ GUI-GENERATED INPUT FILE
C=======================================================================
SUBROUTINE READ_GUI_INPUT(INNAM,ERRNAM,Z1,M1,E0,ESIG,ALPHA,
& ALPHASIG,EF,ESB,SHEATH,ERC,NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,
& KK0R,KDEE1,KDEE2,IPOT,IPOTR,IRL,NLayers,CHEM,DX,RHO,CK,NJ,ZT,
& MT,CO,SBE,ED,BE,CH1,CH2,CH3,CH4,CH5)
IMPLICIT NONE
INTEGER MAXNL,MAXEL
PARAMETER (MAXNL=60)
PARAMETER (MAXEL=12)
CHARACTER*(*) INNAM,ERRNAM
CHARACTER*50 CHEM(MAXNL)
INTEGER I,J,NLAYERS,NH
INTEGER*4 KK0,KK0R,KDEE1,KDEE2,IPOT,IPOTR,IRL
INTEGER*4 NJ(MAXNL)
REAL*8 Z1,M1,E0,ESIG,ALPHA,ALPHASIG
REAL*8 EF,ESB,SHEATH,ERC
REAL*8 RI,RI2,RI3,X0,RD,CW,CA
REAL*8 DX(MAXNL),RHO(MAXNL),CK(MAXNL)
REAL*8 ZT(MAXNL,MAXEL),MT(MAXNL,MAXEL)
REAL*8 CO(MAXNL,MAXEL),SBE(MAXNL,MAXEL)
REAL*8 ED(MAXNL,MAXEL),BE(MAXNL,MAXEL)
REAL*8 CH1(MAXNL,MAXEL),CH2(MAXNL,MAXEL)
REAL*8 CH3(MAXNL,MAXEL),CH4(MAXNL,MAXEL)
REAL*8 CH5(MAXNL,MAXEL)
OPEN(UNIT=99,file=errnam,STATUS='replace')
OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359)
C First line: properties of projectile
C Ordered as: Z, mass number (amu), imp. energy, dist. of imp. energy, angle, dist. angles,
READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC
C Second line: simulation related parameters
C Ordered as: Number of particles, seed, seed, seed, initial depth, RD, depth increment, CA, KK0, KDEE1,KDEE2,IPOT
READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT
& ,IPOTR,IRL
C Third line: Number of layers
read(11,66) NLayers
66 format(9x,I4)
IF(NLayers.LT.1.OR.NLayers.GT.MAXNL) THEN
WRITE(99,*) 'ERROR: number of layers out of range: ',NLayers
WRITE(99,*) 'Valid range is 1..',MAXNL
WRITE(*,*) 'ERROR: number of layers out of range: ',NLayers
STOP
ENDIF
C Here we read the NLayers structure
DO I=1,NLayers
C Chemical formula (chem), Thickness (DX), density (RHO), and correction factor
C (CK, it is always 1.0??), number of elements (nel)
C Possibly add the number of elements in the layer
READ(11,*) chem(I),DX(I),RHO(I),CK(I),NJ(I)
IF(NJ(I).LT.1.OR.NJ(I).GT.MAXEL) THEN
WRITE(99,*) 'ERROR: elements in layer out of range'
WRITE(99,*) 'Layer, NJ, MAXEL = ',I,NJ(I),MAXEL
WRITE(*,*) 'ERROR: elements in layer out of range'
WRITE(*,*) 'Layer, NJ, MAXEL = ',I,NJ(I),MAXEL
STOP
ENDIF
C Replace the following lines with read(11,*) (a(i,j), j = 1, 5)
C Atomic numbers
read(11,*) (ZT(I,j), j = 1, NJ(I))
C Mass numbers (amu)
read(11,*) (MT(I,j), j = 1, NJ(I))
C Concentration (stoichiometry)
read(11,*) (CO(I,j), j = 1, NJ(I))
C Surface binding energy
read(11,*) (SBE(I,j), j = 1, NJ(I))
C Displacement energy, always 30 eV??
read(11,*) (ED(I,j), j = 1, NJ(I))
C Bulk binding energy, usually zero
read(11,*) (BE(I,j), j = 1, NJ(I))
C value A-1 of the ziegler tables
read(11,*) (CH1(I,j), j = 1, NJ(I))
C value A-2 of the ziegler tables
read(11,*) (CH2(I,j), j = 1, NJ(I))
C value A-3 of the ziegler tables
read(11,*) (CH3(I,j), j = 1, NJ(I))
C value A-4 of the ziegler tables
read(11,*) (CH4(I,j), j = 1, NJ(I))
C value A-5 of the ziegler tables
read(11,*) (CH5(I,j), j = 1, NJ(I))
ENDDO
1359 CLOSE(UNIT=11)
RETURN
END
C=======================================================================
C UTILITY SUBROUTINES AND FUNCTIONS
C=======================================================================
SUBROUTINE MOMENTS(FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM,
# X1S,X2S,X3S,X4S,X5S,X6S,Y)
IMPLICIT NONE
REAL*8 FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM,
# X1S,X2S,X3S,X4S,X5S,X6S,Y
REAL*8 U,U2,U3,U4,SIGMA3
REAL*8 X3SP,X4SP,X5SP,X6SP
LOGICAL EQUAL
IF(EQUAL(Y,0.D0))GOTO 10
IF(EQUAL(Y,1.D0))GOTO 10
FIM0=X1S/Y
SEM=X2S/Y-FIM0*FIM0
SIGMA=DSQRT(SEM)
SIGMA3=SEM*SIGMA
U=FIM0/SIGMA
U2=U*U
U3=U2*U
U4=U3*U
X3SP=X3S/(Y*SIGMA3)
X4SP=X4S/(Y*SEM*SEM)
X5SP=X5S/(Y*SEM*SIGMA3)
X6SP=X6S/(Y*SIGMA3*SIGMA3)
THM=X3SP-U*(3.D0+U2)
FOM=X4SP-4.D0*U*X3SP+3.D0*U2*(2.D0+U2)
FIM=X5SP-5.D0*U*X4SP+10.D0*U2*X3SP-2.D0*U3*(5.D0+3.D0*U2)
SIM=X6SP-6.D0*U*X5SP+15.D0*U2*X4SP-20.D0*U3*X3SP+ 5.D0*U4*(3.D0+2
& .D0*U2)
DFIM0=SIGMA/DSQRT(Y)
DSEM=SEM*DSQRT(DMAX1(1.D-20,FOM-1.D0)/(Y))
DTHM=DSQRT(DMAX1(1.D-20, (9.D0+8.75D0*THM*THM+2.25D0*THM*THM*FOM-
& 6.D0*FOM-3.D0*THM*FIM+SIM))/Y)
10 CONTINUE
RETURN
END
SUBROUTINE MOMENTN(FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM,
& X1SY,X2SY,X3SY,X4SY,X5SY,X6SY,X1S,X2S,X3S,X4S,X5S,X6S,Y)
IMPLICIT NONE
REAL*8 FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM, X1SY,X2SY
& ,X3SY,X4SY,X5SY,X6SY,X1S,X2S,X3S,X4S,X5S,X6S,Y
REAL*8 X3SP,X4SP,X5SP,X6SP
REAL*8 U,U2,U3,U4,SIGMA3
LOGICAL EQUAL
IF(EQUAL(Y,0.D0))GOTO 10
IF(EQUAL(Y,1.D0))GOTO 10
X1SY=X1S/Y
X2SY=X2S/Y
X3SY=X3S/Y
X4SY=X4S/Y
X5SY=X5S/Y
X6SY=X6S/Y
FIM0=X1SY
SEM=X2SY-X1SY*X1SY
SIGMA=DSQRT(SEM)
SIGMA3=SEM*SIGMA
U=X1SY/SIGMA
U2=U*U
U3=U2*U
U4=U3*U
X3SP=X3SY/SIGMA3
X4SP=X4SY/(SEM*SEM)
X5SP=X5SY/(SEM*SIGMA3)
X6SP=X6SY/(SIGMA3*SIGMA3)
THM=X3SP-U*(3.D0+U2)
FOM=X4SP-4.D0*U*X3SP+3.D0*U2*(2.D0+U2)
FIM=X5SP-5.D0*U*X4SP+10.D0*U2*X3SP-2.D0*U3*(5.D0+3.D0*U2)
SIM=X6SP-6.D0*U*X5SP+15.D0*U2*X4SP-20.D0*U3*X3SP+ 5.D0*U4*(3.D0+2
& .D0*U2)
DFIM0=SIGMA/DSQRT(Y)
DSEM=SEM*DSQRT(DMAX1(1.D-20,FOM-1.D0)/(Y))
DTHM=DSQRT(DMAX1(1.D-20, (9.D0+8.75D0*THM*THM+2.25D0*THM*THM*FOM-
& 6.D0*FOM-3.D0*THM*FIM+SIM))/Y)
10 CONTINUE
RETURN
END
SUBROUTINE MOMENT(X1SY,X2SY,X3SY,X4SY,X5SY,X6SY ,X1S,X2S,X3S,X4S
& ,X5S,X6S,Y)
IMPLICIT NONE
REAL*8 X1SY,X2SY,X3SY,X4SY,X5SY,X6SY,X1S,X2S,X3S,X4S,X5S,X6S,Y
LOGICAL EQUAL
C Convert accumulated raw sums into per-particle moments by dividing
C through the population size Y. The higher-level MOMENTS/MOMENTN
C routines use these normalized values to build means, spreads and
C uncertainties for the reported observables.
IF(EQUAL(Y,0.D0))GOTO 10
X1SY=X1S/Y
X2SY=X2S/Y
X3SY=X3S/Y
X4SY=X4S/Y
X5SY=X5S/Y
X6SY=X6S/Y
10 RETURN
END
SUBROUTINE DIRCOS(COSX,COSY,COSZ,SINE,CPSI,SPSI,CPHI,SPHI,N)
IMPLICIT NONE
INTEGER N,IV
REAL*8 COSX(N),COSY(N),COSZ(N),SINE(N),CPSI(N),SPSI(N) ,CPHI(N)
& ,SPHI(N)
REAL*8 SRAT,CX2,CY2,CZ2,UNIT
C Rotate the current direction cosines by the sampled scattering
C angles (psi,phi). The output direction is explicitly renormalized
C because repeated collisions amplify small floating-point drift.
DO IV=1,N
SRAT=SPSI(IV)/SINE(IV)
CX2=CPSI(IV)*COSX(IV)+SPSI(IV)*SINE(IV)*CPHI(IV)
CY2=CPSI(IV)*COSY(IV)-SRAT*(COSY(IV)*COSX(IV)*CPHI(IV)-COSZ(IV)
& *SPHI(IV))
CZ2=CPSI(IV)*COSZ(IV)-SRAT*(COSZ(IV)*COSX(IV)*CPHI(IV)+COSY(IV)
& *SPHI(IV))
UNIT = 1.0D0/DSQRT(CX2**2+CY2**2+CZ2**2)
COSX(IV)=CX2*UNIT
COSY(IV)=CY2*UNIT
C MAKE SURE COSZ.NE.0.
COSZ(IV)=DSIGN(DABS(CZ2*UNIT)+1.D-12,CZ2)
SINE(IV)=DSQRT(COSY(IV)*COSY(IV)+COSZ(IV)*COSZ(IV))
ENDDO
RETURN
END
SUBROUTINE VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,N)
IMPLICIT NONE
INTEGER n,I
C
C FETCH A NEW VELOCITY FROM A MAXWELLIAN FLUX AT A SURFACE
C Vectorized version used when multiple projectiles are initialized
C together with thermal surface velocities.
C
REAL*8 FG(2*N),FFG(N),E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
C DIMENSIOM E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
REAL*8 M1,VELC,ZARG
REAL*8 VELX,VELY,VELZ,VELQ,VEL
COMMON/A/ M1,VELC,ZARG
CALL FGAUSS(FG,2*N,N,FFG,N)
DO I=1,N
VELX=DSQRT((FFG(I)*ZARG)**2+VELC)
VELY=FG(I)*ZARG
VELZ=FG(I+N)*ZARG
VELQ=VELX*VELX+VELY*VELY+VELZ*VELZ
VEL=DSQRT(VELQ)
COSX(I)=VELX/VEL
COSY(I)=VELY/VEL
COSZ(I)=VELZ/VEL
SINE(I)=DSQRT(1.D0-COSX(I)*COSX(I))
E(I)=M1*VELQ
ENDDO
RETURN
END
SUBROUTINE VELOC(E,COSX,COSY,COSZ,SINE)
C
C FETCH A NEW VELOCITY FROM A MAXWELLIAN FLUX AT A SURFACE
C Scalar version that reuses Gaussian deviates generated by FGAUSS.
C
IMPLICIT NONE
INTEGER INIV1,INIV3
REAL*8 FG(128),FFG(64)
C BUGFIX: VELOC reuses buffered Gaussian deviates between calls.
C The buffer indices and arrays must be saved and initialized.
SAVE INIV1,INIV3,FG,FFG
DATA INIV1/0/,INIV3/0/
REAL*8 COSX,COSY,COSZ,SINE
REAL*8 M1,VELC,ZARG
REAL*8 VELX,VELY,VELZ,VELQ,VEL,E
COMMON/A/ M1,VELC,ZARG
IF (INIV1.EQ.0) CALL FGAUSS(FG,INIV1,64,FFG,INIV3)
VELX=FFG(INIV3)*ZARG
VELY=FG(INIV1)*ZARG
VELZ=FG(INIV1-1)*ZARG
IF (VELC.GT.0.) THEN
VELX=DSQRT(VELC+VELX**2)
ENDIF
INIV1=INIV1-2
INIV3=INIV3-1
VELQ=VELX*VELX+VELY*VELY+VELZ*VELZ
VEL=DSQRT(VELQ)
COSX=VELX/VEL
COSY=VELY/VEL
COSZ=VELZ/VEL
SINE=DSQRT(1.D0-COSX*COSX)
E=M1*VELQ
RETURN
END
SUBROUTINE FGAUSS (FG,IND,IANZ,FFG,IND2)
C
C RETURN IANZ PAIRS OF RANDOM NUMBER FROM A GAUSSIAN, I.E. IANZ*2
C NUMBERS, AND RETURN THEM IN THE ARRAY FG
C
C THIS FUNCTION SAMPLES FROM A GAUSSIAN OF THE
C FORM DEXP(-(X-ZA)**2/(2.*ZS**2))/(ZS*DSQRT(2*PI))
C ZA=0.
C ZS=1.
C
C IT IS THE BOX-MUELLER METHOD
C FG returns symmetric normal deviates, while FFG stores the positive
C half-space distribution used for Maxwellian flux sampling.
C
IMPLICIT NONE
INTEGER IND,IND2,IANZ,JJ
INTEGER*4 ISEED
REAL*8 PI2,ZZ,ZSIN,ZCOS,AR,ZT
REAL*8 FG(1),FFG(1)
real*4 random(1), ran2(2)
DATA PI2/6.28318530717598D0/
IND=IANZ+IANZ
DO JJ=1,IANZ
C 1. COMPUTE THE SINE AND COSINE OF 2*PI*RAN(1)
C
call ranlux(ran2, 2)
ZZ=PI2*DBLE(ran2(1))
ZSIN=DSIN(ZZ)
ZCOS=DCOS(ZZ)
AR=DLOG(DBLE(ran2(2)))
ZT=DSQRT(-1.0D0*(AR+AR))
FG(JJ+IANZ)=ZT*ZSIN
FG(JJ)=ZT*ZCOS
ENDDO
C
C RETURN IANZ RANDOM NUMBERS FROM A GAUSSIAN FLUX IN THE ARRAY FFG
C
IND2=IANZ
DO JJ=1,IANZ
call ranlux(random, 1)
AR=DLOG(DBLE(random(1)))
FFG(JJ)=DSQRT(-1.D0*(AR+AR))
ENDDO
RETURN
END
SUBROUTINE ENERGV(FE,E,COSX,COSY,COSZ,SINE,N)
C
C FETCH A NEW ENERGY FROM A MAXWELLIAN FLUX AT A SURFACE
C
IMPLICIT NONE
INTEGER N,I
REAL*8 FE(N),E(N),COSX(N),COSY(N),COSZ(N),SINE(N)
REAL*8 M1,EMT
REAL*8 TI,SHEATH,CALFA
COMMON/B/ TI,SHEATH,CALFA
CALL EMAXW(FE,N)
DO I=1,N
EMT=TI*FE(I)**2
COSX(I) = DSQRT((EMT*CALFA*CALFA +SHEATH)/(EMT +SHEATH))
SINE(I) = DSQRT( 1.D0 -COSX(I)*COSX(I))
COSY(I) = SINE(I)
COSZ(I) = 0.D0
E(I) = EMT + SHEATH
ENDDO
RETURN
END
SUBROUTINE ENERG(E,COSX,COSY,COSZ,SINE)
C
C FETCH A NEW ENERGY FROM A MAXWELLIAN FLUX AT A SURFACE
C
IMPLICIT NONE
REAL*8 FE(16)
REAL*8 M1,COSX,SINE,COSY,COSZ,E,EMT
REAL*8 TI,SHEATH,CALFA
COMMON/B/ TI,SHEATH,CALFA
CALL EMAXW(FE,16)
EMT=TI*FE(9)**2
COSX = DSQRT((EMT*CALFA*CALFA +SHEATH)/(EMT +SHEATH))
SINE = DSQRT( 1.D0 -COSX*COSX)
COSY = SINE
COSZ = 0.D0
E = EMT + SHEATH
RETURN
END
SUBROUTINE EMAXW (FE,NUMB)
C
C THIS FUNCTION SAMPLES FROM A MAXWELLIAN (ENERGY) OF THE
C FORM X**2*DEXP(-X**2)*4/DSQRT(PI))
C
C MONTE CARLO SAMPLER C29 (EVERETT, CASHWELL)
C
implicit none
INTEGER NUMB,I
INTEGER*4 ISEED
REAL*8 FE(1)
REAL*8 PI,AR1,AR2
real*4 random(3)
DATA PI/3.14159265358979D0/
DO I=1,NUMB
call ranlux(random, 3)
AR1=DLOG(DBLE(random(1)))
AR2=DLOG(DBLE(random(2)))*(DCOS(PI*0.5*DBLE(random(3))))**2
FE(I)=DSQRT(-1.D0*(AR1+AR2))
ENDDO
RETURN
END
REAL*8 FUNCTION CVMGT(A, B, C)
C this function returns A if C is tuer and B otherwise
C why is it needed??
IMPLICIT NONE
REAL*8 A,B
LOGICAL C
CVMGT = B
IF ( C ) CVMGT = A
RETURN
END
SUBROUTINE SCOPY(IM,A,INCA,B,INCB)
IMPLICIT NONE
INTEGER*4 INCA,INCB,IM,JA,JB
INTEGER J
REAL*8 A(*),B(*)
JA = IM * IABS(INCA)
IF (INCA .GT. 0) JA = 1
JB = IM * IABS(INCB)
IF (INCB .GT. 0) JB = 1
DO J = 1,IM
B(JB) = A(JA)
JA = JA + INCA
JB = JB + INCB
ENDDO
RETURN
END
FUNCTION ILLZ(N,A,K)
IMPLICIT NONE
LOGICAL A(*)
INTEGER K,L,N,I
INTEGER*4 ILLZ
IF(K.GT.0) THEN
L=N+1
DO I=N,1,-1
IF(A(I)) L=I
ENDDO
ELSE
L=0
DO I=1,N
IF(A(I)) L=I
ENDDO
L=N+1-L
ENDIF
ILLZ=L-1
RETURN
END
INTEGER FUNCTION ISRCHFGE(N,ARRAY,INC,TARGT)
IMPLICIT NONE
INTEGER I,N,J,INC
REAL*8 ARRAY(N),TARGT
J=1
IF(INC.LT.0) J=N*(-INC)
DO I=1,N
IF(ARRAY(J).GE.TARGT) GO TO 200
J=J+INC
ENDDO
200 ISRCHFGE=I
RETURN
END
INTEGER FUNCTION ISRCHFGT(N,ARRAY,INC,TARGT)
IMPLICIT NONE
INTEGER I,N,J,INC
REAL*8 ARRAY(N),TARGT
C Return the first index whose value is strictly greater than TARGT.
C This is used to map a continuous coordinate/depth onto the current
C layer or histogram bin in monotonic arrays.
J=1
IF(INC.LT.0) J=N*(-INC)
DO I=1,N
IF(ARRAY(J).GT.TARGT) GO TO 200
J=J+INC
ENDDO
200 ISRCHFGT=I
RETURN
END
SUBROUTINE ENERGGAUSS(ISEED2,Esig,Epar,E0)
IMPLICIT NONE
INTEGER*4 ISEED2
REAL*8 E0,Esig,Epar
REAL*8 p1,p2,p3,pi
real*4 random(2)
DATA pi/3.14159265358979D0/
C Sample a Gaussian-distributed incident energy around E0 using a
C Box-Mueller transform. Epar is the per-projectile starting energy.
call ranlux(random, 2)
p1 = Esig*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1))))
p2 = 2.D0*pi*DBLE(random(2))
p3 = p1*DCOS(p2)
Epar= E0-p3
RETURN
END
SUBROUTINE ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, CALFA
& ,SALFA,BW)
IMPLICIT NONE
INTEGER*4 ISEED3
REAL*8 ALPHA,ALPHASIG,ALPHApar
REAL*8 BW,ALFA,CALFA,SALFA
REAL*8 p1,p2,p3,pi
real*4 random(2)
DATA pi/3.14159265358979D0/
C Sample a Gaussian-distributed incidence angle around ALPHA and
C return both the angle in radians (ALFA) and its sine/cosine.
call ranlux(random, 2)
p1 = ALPHASIG*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1))))
p2 = 2.D0*pi*DBLE(random(2))
p3 = p1*DCOS(p2)
ALPHApar= ALPHA-p3
IF(ALPHApar.LT.0.D0) THEN
ALPHApar=DABS(ALPHApar)
ENDIF
ALFA = ALPHApar/BW
CALFA = DCOS(ALFA)
SALFA = DSIN(ALFA)
RETURN
END
LOGICAL FUNCTION EQUAL(F1,F2)
IMPLICIT NONE
REAL*8 F1,F2
REAL*8 TINY
PARAMETER (TINY = 1.0D-10)
C Centralized floating-point comparison used throughout the transport
C code to avoid exact equality checks on derived real-valued quantities.
IF (DABS(F1-F2).LE.TINY) THEN
EQUAL = .TRUE.
ELSE
EQUAL = .FALSE.
ENDIF
RETURN
END
SUBROUTINE TimeStamp(day,month,year,hour,min,sec,seconds_total)
C Subroutine to return date and time for start/stop of simulation
IMPLICIT NONE
INTEGER Date_Time(8)
INTEGER*4 days_total
INTEGER*4 seconds_total
CHARACTER Real_Clock(3)*12
CHARACTER month*4,day*2
CHARACTER year*4,hour*2
CHARACTER min*2,sec*2
CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3),
& Date_Time)
IF(Date_Time(2).EQ.1) THEN
month='Jan.'
days_total=Date_Time(3)
ELSEIF(Date_Time(2).EQ.2) THEN
month='Feb.'
days_total=31+Date_Time(3)
ELSEIF(Date_Time(2).EQ.3) THEN
month='Mar.'
days_total=31+28+Date_Time(3)
ELSEIF(Date_Time(2).EQ.4) THEN
month='Apr.'
days_total=31+28+31+Date_Time(3)
ELSEIF(Date_Time(2).EQ.5) THEN
month='May '
days_total=31+28+31+30+Date_Time(3)
ELSEIF(Date_Time(2).EQ.6) THEN
month='Jun.'
days_total=31+28+31+30+31+Date_Time(3)
ELSEIF(Date_Time(2).EQ.7) THEN
month='Jul.'
days_total=31+28+31+30+31+30+Date_Time(3)
ELSEIF(Date_Time(2).EQ.8) THEN
month='Aug.'
days_total=31+28+31+30+31+30+31+Date_Time(3)
ELSEIF(Date_Time(2).EQ.9) THEN
month='Sep.'
days_total=31+28+31+30+31+30+31+31+Date_Time(3)
ELSEIF(Date_Time(2).EQ.10) THEN
month='Oct.'
days_total=31+28+31+30+31+30+31+31+30+Date_Time(3)
ELSEIF(Date_Time(2).EQ.11) THEN
month='Nov.'
days_total=31+28+31+30+31+30+31+31+30+31+Date_Time(3)
ELSE
month='Dec.'
days_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3)
ENDIF
C in seconds from beginning of year
seconds_total=Date_Time(7)+(Date_Time(6)*60)+(Date_Time(5) *3600)
& +(days_total-1)*86400
READ(Real_Clock(1)(1:4),'(A4)')year
READ(Real_Clock(1)(7:8),'(A2)')day
READ(Real_Clock(2)(1:2),'(A2)')hour
READ(Real_Clock(2)(3:4),'(A2)')min
READ(Real_Clock(2)(5:6),'(A2)')sec
RETURN
END
C=======================================================================
C RANDOM NUMBER GENERATOR
C=======================================================================
SUBROUTINE RANLUX(RVEC,LENV)
C Subtract-and-borrow random number generator proposed by
C Marsaglia and Zaman, implemented by F. James with the name
C RCARRY in 1991, and later improved by Martin Luescher
C in 1993 to produce "Luxury Pseudorandom Numbers".
C Fortran 77 coded by F. James, 1993
C
C LUXURY LEVELS.
C ------ ------ The available luxury levels are:
C
C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
C and Zaman, very long period, but fails many tests.
C level 1 (p=48): considerable improvement in quality over level 0,
C now passes the gap test, but still fails spectral test.
C level 2 (p=97): passes all known tests, but theoretically still
C defective.
C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
C correlations have very small chance of being observed.
C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
C
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C!!! Calling sequences for RANLUX: ++
C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
C!!! 32-bit random floating point numbers between ++
C!!! zero (not included) and one (also not incl.). ++
C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
C!!! one 32-bit integer INT and sets Luxury Level LUX ++
C!!! which is integer between zero and MAXLEV, or if ++
C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
C!!! should be set to zero unless restarting at a break++
C!!! point given by output of RLUXAT (see RLUXAT). ++
C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
C!!! which can be used to restart the RANLUX generator ++
C!!! at the current point by calling RLUXGO. K1 and K2++
C!!! specify how many numbers were generated since the ++
C!!! initialization with LUX and INT. The restarting ++
C!!! skips over K1+K2*E9 numbers, so it can be long.++
C!!! A more efficient but less convenient way of restarting is by: ++
C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
C!!! 32-bit integer seeds, to be used for restarting ++
C!!! ISVEC must be dimensioned 25 in the calling program ++
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
DIMENSION RVEC(LENV)
DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
PARAMETER (MAXLEV=4, LXDFLT=3)
DIMENSION NDSKIP(0:MAXLEV)
DIMENSION NEXT(24)
PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
PARAMETER (ITWO24=2**24, ICONS=2147483563)
SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
INTEGER LUXLEV
LOGICAL NOTYET
DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
DATA I24,J24,CARRY/24,10,0./
C default
C Luxury Level 0 1 2 *3* 4
DATA NDSKIP/0, 24, 73, 199, 365 /
Corresponds to p=24 48 97 223 389
C time factor 1 2 3 6 10 on slow workstation
C 1 1.5 2 3 5 on fast mainframe
C
C NOTYET is .TRUE. if no initialization has been performed yet.
C Default Initialization by Multiplicative Congruential
IF (NOTYET) THEN
NOTYET = .FALSE.
JSEED = JSDFLT
INSEED = JSEED
LUXLEV = LXDFLT
NSKIP = NDSKIP(LUXLEV)
LP = NSKIP + 24
IN24 = 0
KOUNT = 0
MKOUNT = 0
TWOM24 = 1.
DO 25 I= 1, 24
TWOM24 = TWOM24 * 0.5
K = JSEED/53668
JSEED = 40014*(JSEED-K*53668) -K*12211
IF (JSEED .LT. 0) JSEED = JSEED+ICONS
ISEEDS(I) = MOD(JSEED,ITWO24)
25 CONTINUE
TWOM12 = TWOM24 * 4096.
DO 50 I= 1,24
SEEDS(I) = REAL(ISEEDS(I))*TWOM24
NEXT(I) = I-1
50 CONTINUE
NEXT(1) = 24
I24 = 24
J24 = 10
CARRY = 0.
IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
ENDIF
C
C The Generator proper: "Subtract-with-borrow",
C as proposed by Marsaglia and Zaman,
C Florida State University, March, 1989
C
DO 100 IVEC= 1, LENV
UNI = SEEDS(J24) - SEEDS(I24) - CARRY
IF (UNI .LT. 0.) THEN
UNI = UNI + 1.0
CARRY = TWOM24
ELSE
CARRY = 0.
ENDIF
SEEDS(I24) = UNI
I24 = NEXT(I24)
J24 = NEXT(J24)
RVEC(IVEC) = UNI
C small numbers (with less than 12 "significant" bits) are "padded".
IF (UNI .LT. TWOM12) THEN
RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
C and zero is forbidden in case someone takes a logarithm
IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
ENDIF
C Skipping to luxury. As proposed by Martin Luscher.
IN24 = IN24 + 1
IF (IN24 .EQ. 24) THEN
IN24 = 0
KOUNT = KOUNT + NSKIP
DO 90 ISK= 1, NSKIP
UNI = SEEDS(J24) - SEEDS(I24) - CARRY
IF (UNI .LT. 0.) THEN
UNI = UNI + 1.0
CARRY = TWOM24
ELSE
CARRY = 0.
ENDIF
SEEDS(I24) = UNI
I24 = NEXT(I24)
J24 = NEXT(J24)
90 CONTINUE
ENDIF
100 CONTINUE
KOUNT = KOUNT + LENV
IF (KOUNT .GE. IGIGA) THEN
MKOUNT = MKOUNT + 1
KOUNT = KOUNT - IGIGA
ENDIF
RETURN
C
C Entry to input and float integer seeds from previous run
ENTRY RLUXIN(ISDEXT)
NOTYET = .FALSE.
TWOM24 = 1.
DO I= 1, 24
NEXT(I) = I-1
TWOM24 = TWOM24 * 0.5
ENDDO
NEXT(1) = 24
TWOM12 = TWOM24 * 4096.
WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
WRITE(6,'(5X,5I12)') ISDEXT
DO 200 I= 1, 24
SEEDS(I) = REAL(ISDEXT(I))*TWOM24
200 CONTINUE
CARRY = 0.
IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
ISD = IABS(ISDEXT(25))
I24 = MOD(ISD,100)
ISD = ISD/100
J24 = MOD(ISD,100)
ISD = ISD/100
IN24 = MOD(ISD,100)
ISD = ISD/100
LUXLEV = ISD
IF (LUXLEV .LE. MAXLEV) THEN
NSKIP = NDSKIP(LUXLEV)
WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
+ LUXLEV
ELSE IF (LUXLEV .GE. 24) THEN
NSKIP = LUXLEV - 24
WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
ELSE
NSKIP = NDSKIP(MAXLEV)
WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
LUXLEV = MAXLEV
ENDIF
INSEED = -1
RETURN
C
C Entry to ouput seeds as integers
ENTRY RLUXUT(ISDEXT)
DO 300 I= 1, 24
ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
300 CONTINUE
ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
RETURN
C
C Entry to output the "convenient" restart point
ENTRY RLUXAT(LOUT,INOUT,K1,K2)
LOUT = LUXLEV
INOUT = INSEED
K1 = KOUNT
K2 = MKOUNT
RETURN
C
C Entry to initialize from one or three integers
ENTRY RLUXGO(LUX,INS,K1,K2)
IF (LUX .LT. 0) THEN
LUXLEV = LXDFLT
ELSE IF (LUX .LE. MAXLEV) THEN
LUXLEV = LUX
ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
LUXLEV = MAXLEV
WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
ELSE
LUXLEV = LUX
DO 310 ILX= 0, MAXLEV
IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
310 CONTINUE
ENDIF
IF (LUXLEV .LE. MAXLEV) THEN
NSKIP = NDSKIP(LUXLEV)
WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
+ LUXLEV,' P=', NSKIP+24
ELSE
NSKIP = LUXLEV - 24
WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
ENDIF
IN24 = 0
IF (INS .LT. 0) WRITE (6,'(A)')
+ ' Illegal initialization by RLUXGO, negative input seed'
IF (INS .GT. 0) THEN
JSEED = INS
WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
+ JSEED, K1,K2
ELSE
JSEED = JSDFLT
WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
ENDIF
INSEED = JSEED
NOTYET = .FALSE.
TWOM24 = 1.
DO 325 I= 1, 24
TWOM24 = TWOM24 * 0.5
K = JSEED/53668
JSEED = 40014*(JSEED-K*53668) -K*12211
IF (JSEED .LT. 0) JSEED = JSEED+ICONS
ISEEDS(I) = MOD(JSEED,ITWO24)
325 CONTINUE
TWOM12 = TWOM24 * 4096.
DO 350 I= 1,24
SEEDS(I) = REAL(ISEEDS(I))*TWOM24
NEXT(I) = I-1
350 CONTINUE
NEXT(1) = 24
I24 = 24
J24 = 10
CARRY = 0.
IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
C If restarting at a break point, skip K1 + IGIGA*K2
C Note that this is the number of numbers delivered to
C the user PLUS the number skipped (if luxury .GT. 0).
KOUNT = K1
MKOUNT = K2
IF (K1+K2 .NE. 0) THEN
DO 500 IOUTER= 1, K2+1
INNER = IGIGA
IF (IOUTER .EQ. K2+1) INNER = K1
DO 450 ISK= 1, INNER
UNI = SEEDS(J24) - SEEDS(I24) - CARRY
IF (UNI .LT. 0.) THEN
UNI = UNI + 1.0
CARRY = TWOM24
ELSE
CARRY = 0.
ENDIF
SEEDS(I24) = UNI
I24 = NEXT(I24)
J24 = NEXT(J24)
450 CONTINUE
500 CONTINUE
C Get the right value of IN24 by direct calculation
IN24 = MOD(KOUNT, NSKIP+24)
IF (MKOUNT .GT. 0) THEN
IZIP = MOD(IGIGA, NSKIP+24)
IZIP2 = MKOUNT*IZIP + IN24
IN24 = MOD(IZIP2, NSKIP+24)
ENDIF
C Now IN24 had better be between zero and 23 inclusive
IF (IN24 .GT. 23) THEN
WRITE (6,'(A/A,3I11,A,I5)')
+ ' Error in RESTARTING with RLUXGO:',' The values', INS,
+ K1, K2, ' cannot occur at luxury level', LUXLEV
IN24 = 0
ENDIF
ENDIF
RETURN
END
Reference in New Issue View Git Blame Copy Permalink