diff --git a/trimsp/src/trimspNL.F b/trimsp/src/trimspNL.F new file mode 100644 index 0000000..716c0d2 --- /dev/null +++ b/trimsp/src/trimspNL.F @@ -0,0 +1,4855 @@ +C Version TrimSpNL ----> N Layer +C +C Created Juni 2000 ---- Testversion +C +C *** +C * C * COPYRIGHT W.ECKSTEIN, IPP GARCHING, FRG +C *** +C +C +C +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 MAJOR CHANGES: +C +C Sep 1998: to create an executable for PC running under WIN95 +C using the DIGITAL VISUAL FORTRAN compiler (ed. Aug. 97) +C Insertion of !DEC$REAL:8 (all REAL are REAL*8) +C Conversion of function to double precession +C Conversion of function to integer*4 +C Conversion of random number generator from DRAND48() +C to DBLE(RAN(ISEED)) +C UNIT 17 (data to tape) disabled +C UNIT 11 (input data) changed to file name EINGABE1.inp +C UNIT 21 (output.data) changed to file name AUSGABE1.out +C Insertion of UNIT 22 (only range data), file name AUSGABE1.rge +C Dec 1998: Introduction of gaussian distributed projectile energies +C new variable RI2 = random number initializer for gaussian energy distribution +C new variable ISEED2 = random number for gaussian energy distribution +C new variable Esig = sigma of the energy distribution +C if Esig=0. then fixed projectile energy +C new variable Epar = projectile energy +C new subroutine ENERGGAUSS = for the gaussian distribution +C new variables in subroutine p1,p2,p3, necessary for calculation of +C the gaussian energy distribution +C Mar 1999: LOGICAL FUNCTION EQUAL inserted. This function is used for comparision +C of REAL*8 variables with 0. +C OUTPUT to files AUSGABE1: +C bug fixed for the output of the CH(i), bug was also present +C in the original and the older version of TRIMSPP*C +C if OUTPUT file exists then program asks for new OUTPUT filename +C all variables are now defined as REAL*8 or INTEGER*4 +C for variable conversion from REAL to INTEGER or INTEGER to REAL +C conversion function inserted (IDINT, DFLOTJ) +C use of correct MIN or MAX functions +C data initialization for all variables (like in BHam) but without +C DATA ier/102*0/ give an error in calculation of ions stopped in layer 2 - 3 +C Unit 99 for file ausgabe1.err inserted, if and IF... (see below is true +C then a message is included into this file +C IF's inserted - can be found by using find WRITE(99, +C DABS inserted for DLOG and DLOG10 arguments +C Mai 1999: Version h +C for TRIMSP simulations running in batch mode (i.e. to calculate +C a serie of different range profiles using the same target an output +C file FOR33 is created. In this file the following parameters are +C written: +C E0,Esig,NH,IIM,IB,IT,IRPL(1),IRPL(2),IRPL(3),fix0,sigmax +C the FILE open statement is inserted after label 1502 +C new variable column(100)*104 and colcount inserted. +C Mai 1999: Version i +C tried to fix bug for variable C2 (in Do 65 loop). +C Sometimes this variable becames smaller then 10-10 (an underflow occurs, +C the variable C2 is NaN, variables S2 and DENS, which are directly calculated +C from C2 are also NaN. If variable C2 is less then 10-10 then C2 is 10-10. +C S2 and DENS are calculated with the resetted value of C2. +C Range file (output to unit 22) slightly changed. Now the midth of a channel +C is calculated and only the number of particles is written to the file +C +C +C Jun 1999: Version j1 +C changed calculation of Firsoc screening length according to the IRCU 49 report. +C the Andersen-Ziegler tables (file stopping.dat) are revised,i.e. the coefficients +C from the ICRU 49 report are included for all elements (file stopicru.dat) +C +C PROGRAM DATMAK2c: this program creates the input files for TRIMSPP4*. +C with this program one can change the energy, +C the energy distribution, the number of muons, and the +C layer thickness easily. This is for making different +C simulations (i.e. energy scan, sigma scan, particle +C number scan, layer thickness scan). A corresponding +C batch file is created with this program too. +C Nov 1999: Version j1a +C UNIT33: now Etrans, sigmaEtrans, Eback, sigmaEback included +C Bug fixed in line 2084, location of 'write to unit33' changed. +C Dec 1999: Version k +C inclusion of variables +C tryE :how often a random energy is calculated by the random generator +C for large E0 this number should be the same as the number of projectiles nh +C negE :how often a negative energy is calculated by the random generator +C tryE - negE should be nh +C both variables are type integer, both variables included in the output file +C and in the fort33.file +C inclusion of alphasig: one dimensional distribution of the angle of incidence +C if alphasig .NE. 0 THEN a gaussian distribution of alpha will be calculated +C in the subroutine ALPHAGAUSS +C if the calculated angle is < 0 then the absolute value is taken because of +C the options alpha = -1 or alpha = -2 (see file TRVMC95-v4k.txt) +C inclusion of RI3: random seed for the calculation of the gaussian distribution +C of alpha. +C Header line for file for33 included. +C Jan 2000: no new version but +C included energy scaling of particle reflection coefficients after Thomas et al. +C NIM B69 (1992) 427 +C the calculated particle reflection coefficients prcoeff are written to the fort.33 file +C prc is calculated if 1st layer consists only one element ! +C included constants prc(1) - prc(6) +C Thomas-Fermi reduced energy: epsilon +C output of E0 and Esig now in keV +C Apr 2000: BUG fixed, Var CHM2 removed, i.e. stopping power for energies below 10 keV +C was a little bit overestimated. now version TrimSpP4L +C Prg Datmak4L creates input files. +C Jun 2000: Source code enlarged to simulate implantation profiles in up to 7 layers. +C Each layer can have up to 5 different elements. +C BUG fixed in line 1470,i.e. in calculation of scattering angle for the +C Moliere potential the variable rrr1 was not handled correct. The length of +C line 1470 was too long, RRR1 was read as RRR which has the value 0. +C Due to the fact, that all calculation based on the other interaction potentials +C (Krypton Carbon and ZBL) give more or less the same inplantation profiles, this +C bug has a minor influence to the implantation profile. +C Variables DateTime(8) = Integer array and +C Real_Clock(3) = Character array included for date and time output to file +C Changes for output file to UNIT 33: +C Size of character variable column enlarged from 214 to 246 +C Program datmak7a creates input files for this TrimSp release. +C true calculation, i.e. indendent of bin width, of particles stopped +C in layers - see UpTiefe,LowTiefe,number_in_layer(7) +C INTEGER check_layer_flag : IF 1 then calculated implantation profile +C in the 100 depth intervals agrees with +C number of particles stopped in the different layers +C IF 0 not, message will be written in the range file +C UNIT 21 and UNIT22 +c Nov 2000: Tanya Riseman +c Starts porting program to Open VMS and DEC UNIX. +c Compile options: f90/list/warn/ext +c Minor problems with continuation characters in wrong column and +c with a few variables undeclared, mostly functions. +c Make code fit on 72 columns, because -extend_source +c does not appear to work on DEC Unix. +c Start changing strings in format and write statments so that +c strings don't straddle continuation character in +c column 6. Straddling can be non-portable! +c Comment out !DEC$REAL:8 because all +c reals are already declared REAL*8. Add implicit none. +c +c Jun 2002: Thomas Prokscha PSI +c Stopped porting to f90, use f77 only. +c replaced all non-standard f77 function by standard functions. +c use ranlux random number generator from the CERN library (libmathlib.a must +c be installed) to get rid of machine specific random number generators. +c Add pre-compiler instructions for making different output for the .rge files +c on Windows and Unix. +c +c Oct 2002: Thomas Prokscha PSI +c corrected error in the calculation of the Thomas-Fermi reduced energy: +c it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3)) +c it must be: +c 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3))) +c +c Aug 2009: Zaher Salman PSI +c Changed input file reading to non-formatted. This still works +c with old input files, but makes it much easier to create new +c input files without the hassle of formatting. +c Included the source of ranlux for random numbers generation +c into trimsp7l source. No need for cern libraries to be installed. +c +c Feb 2013: Zaher Salman PSI +c Started cleaning up the code. +c When possible remove loop numbres and use do-enddo instead. +c Using proper fortran line indentation (emacs style). +c Created TimeStamp subroutine to generate time stamp and removed +c original code from main (two places) +c +c------------------------------------------- +c check OS +c +#if defined( _WIN32 ) +#define OS_WIN +#else +#define OS_UNIX +#endif + + IMPLICIT NONE + LOGICAL TEST(64),TESTR(2000),TEST1(2000) + LOGICAL EQUAL + INTEGER*4 ISRCHEQ,ISRCHFGT,ISRCHFGE,ILLZ + INTEGER N,L,LL,NH,NUM,KK + INTEGER I,J,IV + INTEGER tryE,negE + INTEGER COLCOUNT + INTEGER depth_interval_flag + INTEGER OldNew + 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(7),laufzahl + INTEGER*4 IRP(0:101),IPL(100),IPLB(100),IPLT(100) + INTEGER*4 ICD(100,35),ICDT(100),ICDJT(35),ICDIRJ(35,35),ICDR(100 + & ,35),ICDTR(100),ICDJTR(35) ,ICDIRI(100,35,35),ICDIRN(100,35) + & ,ICDITR(35) + INTEGER*4 KADB(20),KADT(20),KADS(20),KADST(20) ,KADRIP(20,30) + & ,KADRIS(20,30),KADROP(20,30),KADROS(20,30) ,KADSJ(20,30) + & ,KADSL(20,6),KDSTJ(20,30),KDSTL(20,6) + INTEGER*4 IBSP(35),ISPAL(7),IBSPL(35) ,ISPIP(35),ISPIS(35) + & ,ISPOP(35),ISPOS(35) + INTEGER*4 ITSP(35),ISPALT(7) + & ,ISPIPT(35),ISPIST(35),ISPOPT(35),ISPOST(35) + INTEGER*4 KO(600,35,2) + INTEGER*4 MEAB(102,22),MAGB(62,22),MEAGB(102,36,22) ,MEABL(75,21) + & ,MEPB(102,102) + INTEGER*4 MEAT(102,22),MAGT(62,22),MEAGT(102,36,22), + & MEATL(75,21),MEPT(102,102) + INTEGER*4 MEAS(102,22,30),MAGS(62,22,30),MAGSA(62,32,30) + & ,MEAGS(102,12,22,30) ,MEASL(75,21,30) + INTEGER*4 MEAST(102,22,30),MAGST(62,22,30) ,MEASTL(75,21,30) + INTEGER*4 NJ(7),JT(7),ILD(7) + INTEGER*4 LLL(64),JJJ(64),IK(64) + INTEGER*4 me(5000),nli(0:7),irpl(7) + 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 LMAX,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:101) ,CASMOT(100),PHON(100),DENT(100),ION(100) + & ,DMGN(100) ,CASMOTR(100),PHONR(100),DENTR(100),IONR(100) + & ,DMGNR(100) ,ELGD(100),ELGDR(100) + REAL*8 ELE(100,35),ELI(100,35),ELP(100,35),ELD(100,35) ,ELET(35) + & ,ELIT(35),ELPT(35),ELDT(35) ,ELER(100,35),ELIR(100,35) + & ,ELPR(100,35),ELDR(100,35) ,ELETR(35),ELITR(35),ELPTR(35) + & ,ELDTR(35) + REAL*8 AI(20),RKADB(20),RKADT(20) ,RKADS(20),RKADST(20) + & ,RKADSJ(20,30),RKADSL(20,7) ,RKDSTJ(20,30),RKDSTL(20,7) + REAL*8 EBSP(35),ESPAL(7) ,SPY(35),SPE(35),REY(35),EMSP(35) + & ,ESPIP(35),ESPIS(35),ESPOP(35),ESPOS(35) ,RIP(35),RIS(35) + & ,ROP(35),ROS(35) ,REIP(35),REIS(35),REOP(35),REOS(35) + & ,ESPMIP(35),ESPMIS(35),ESPMOP(35),ESPMOS(35) ,RIPJ(35) + & ,RISJ(35),ROPJ(35),ROSJ(35) ,REIPJ(35),REISJ(35),REOPJ(35) + & ,REOSJ(35) + REAL*8 ETSP(35),ESPALT(7) ,SPYT(35),SPET(35),REYT(35),EMSPT(35) + & ,ESPIPT(35),ESPIST(35),ESPOPT(35),ESPOST(35) ,RIPT(35) + & ,RIST(35),ROPT(35),ROST(35) ,REIPT(35),REIST(35),REOPT(35) + & ,REOST(35) ,ESPMIPT(35),ESPMIST(35),ESPMOPT(35),ESPMOST(35) + REAL*8 SPEM(35),SPE2S(35),SPE3S(35),SPE4S(35),SPE5S(35) ,SPE6S(35) + & ,VSPE(35),SSPE(35),GSPE(35),BSPE(35) + REAL*8 SPE1SL(35),SPE2SL(35),SPE3SL(35),SPE4SL(35),SPE5SL(35) + & ,SPE6SL(35) + REAL*8 ELOG(75),EMA(62,22),EABL(75) + REAL*8 EMAT(62,22),EATL(75),EASL(75,30),EASTL(75,30) + REAL*8 FG(128),FFG(64) + REAL*8 XX(7),DX(7),RHO(7),Z2(7),M2(7),LM(7),PDMAX(7),ARHO(7),AM(7) + & ,FM(7),EPS0(7),ASIG(7),K2(7),CK(7) ,KLM1(7),SB(7),DLI(7) + REAL*8 UpTiefe,LowTiefe + REAL*8 ZT(7,5),MT(7,5),CO(7,5),SBE(7,5),ED(7,5),BE(7,5) ,COM(5,7) + REAL*8 MU(35,35),EC(35,35),A(35,35),F(35,35) ,KL(35,35),KOR(35,35) + & ,KLM(7,35) + REAL*8 MU1(35),EC1(35),A1(35),F1(35),KL1(35),KOR1(35) ,DI(35) + & ,EP(35),ZZ(35),TM(35) + REAL*8 CH1(7,5),CH2(7,5),CH3(7,5),CH4(7,5),CH5(7,5) + REAL*8 CHM1(7) + REAL*8 SM(64),SH(64,5) + REAL*8 FIESB(30),SEESB(30),THESB(30),FOESB(30) ,SGMESB(30) + & ,DFIESB(30),DSEESB(30),DTHESB(30) + 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 fuer part. reflec. coeff. Berechnung + 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 + 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,ran2(2) +C CHARACTER Variables + CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2 + CHARACTER filein*8,inext*4,fileout*8,outext*4,innam*12,outnam*12 + CHARACTER rgenam*12,rgeext*4,errnam*12,errext*4 + CHARACTER errcom*72 + CHARACTER COLUMN(100)*246 + 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 + + COMMON /A/ M1,VELC,ZARG + COMMON /B/ TI,SHEATH,CALFA + + 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'/ + DATA filein/'eingabe1'/,fileout/'ausgabe1'/ + + DATA ET/0.D0/,PLST/0.D0/,PL2ST/0.D0/,PL3ST/0.D0/ + DATA PL4ST/0.D0/,PL5ST/0.D0/,PL6ST/0.D0/ + DATA SEM/0.D0/,IT/0/,NPROJ/0/,NREC1/0/,NREC2/0/ + DATA NH/0/,IB/0/,IBL/0/,NJ/7*0/,NLI/8*0/,DLI/7*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/2244*0/,MEABL/1575*0/,MAGB/1364*0/ + DATA MEPB/10404*0/,MEAGB/80784*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/ + 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/100*0.D0/,DENT/100*0.D0/ + DATA DMGN/100*0.D0/,ION/100*0.D0/,PHON/100*0.D0/ + DATA PHONR/100*0.D0/ + DATA ELGD/100*0.D0/,ELGDR/100*0.D0/ + DATA ICDT/100*0/,ICDTR/100*0/ + DATA ICDR/3500*0/,ICDIRN/3500*0/,IONR/100*0.D0/ + DATA DENTR/100*0.D0/,DMGNR/100*0.D0/ + DATA IPL/100*0/,IPLB/100*0/,IPLT/100*0/ + DATA IRPL/7*0/ + DATA ICDJT/35*0/,ICDJTR/35*0/,ICDITR/35*0/ + DATA ICD/3500*0/,ELP/3500*0.D0/,ELD/3500*0.D0/ + DATA ELE/3500*0.D0/,ELI/3500*0.D0/ + DATA ICDIRI/122500*0/ + DATA ICSUM/0/,ICSUMS/0/,ICDI/0/,ISPA/0/,ISPAT/0/ + DATA Z2/7*0.D0/,M2/7*0.D0/ + DATA KLM1/7*0.D0/,CHM1/7*0.D0/,SB/7*0.D0/,KLM/245*0.D0/ + DATA ME/5000*0/,EMX/0.D0/,ESPAT/0.D0/,ESPA/0.D0/ + DATA IBSP/35*0/,IBSPL/35*0/,EBSP/35*0.D0/,ISPAL/7*0/ + DATA ITSP/35*0/,ETSP/35*0.D0/ + DATA ESPAL/7*0.D0/,ESPALT/7*0.D0/,ISPALT/7*0/ + DATA SPE2S/35*0.D0/,SPE3S/35*0.D0/,SPE4S/35*0.D0/ + DATA SPE5S/35*0.D0/,SPE6S/35*0.D0/ + DATA SPE1SL/35*0.D0/,SPE2SL/35*0.D0/,SPE3SL/35*0.D0/ + DATA SPE4SL/35*0.D0/,SPE5SL/35*0.D0/,SPE6SL/35*0.D0/ + DATA ELET/35*0.D0/,ELPT/35*0.D0/,ELDT/35*0.D0/,ELIT/35*0.D0/ + DATA ELETR/35*0.D0/,ELITR/35*0.D0/,ELPTR/35*0.D0/ + DATA ELDTR/35*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 /7*0/ + + innam=filein//inext + outnam=fileout//outext + rgenam=fileout//rgeext + errnam=fileout//errext + + if (OldNew(innam).eq.1) then + OPEN(UNIT=99,file=errnam) + OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359) + +C This part reads the input file (new format) +C First line: properties of projectile + READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC +C Second line: simulation related parameters + 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(8x,I4) +C Here we read the NLayers structure + DO I=1,NLayers +C Thickness (DX), density (RHO), and correction factor (CK, it is +C always 1.0??) Atomic numbers + READ(11,*) DX(I),RHO(I),CK(I) +C Atomic numbers + READ(11,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) +C Mass numbers (amu) + READ(11,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) +C Concentration + READ(11,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) +C Surface binding energy + READ(11,*) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) +C Displacement energy + READ(11,*) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) +C Bulk binding energy + READ(11,*) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) +C value A-1 of the ziegler tables + READ(11,*) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) +C value A-2 of the ziegler tables + READ(11,*) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) +C value A-3 of the ziegler tables + READ(11,*) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) +C value A-4 of the ziegler tables + READ(11,*) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) +C value A-5 of the ziegler tables + READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) + ENDDO + else + OPEN(UNIT=99,file=errnam) + OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359) + +C This part reads the input file (old format, 7 layers) +C First line: properties of projectile + READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC +C Second line: simulation related parameters + READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT + & ,IPOTR,IRL +C Third line: layer structure. To be replaced by number of layers +C and then each layer with its properties: Thickness (DX), density +C (RHO), and correction factor (CK, it is always 1.0??) + READ(11,*) DX(1),DX(2),DX(3),DX(4),DX(5),DX(6),DX(7),RHO(1) + & ,RHO(2),RHO(3),RHO(4),RHO(5),RHO(6),RHO(7), CK(1),CK(2) + & ,CK(3) ,CK(4),CK(5),CK(6),CK(7) +C Here we read the 7 layer structure + DO I=1,7 +C Atomic numbers + READ(11,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) +C Mass numbers (amu) + READ(11,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) +C Concentration + READ(11,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) +C Surface binding energy + READ(11,*) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) +C Displacement energy + READ(11,*) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) +C Bulk binding energy + READ(11,*) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) +C value A-1 of the ziegler tables + READ(11,*) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) +C value A-2 of the ziegler tables + READ(11,*) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) +C value A-3 of the ziegler tables + READ(11,*) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) +C value A-4 of the ziegler tables + READ(11,*) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) +C value A-5 of the ziegler tables + READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) + ENDDO + endif + + 1359 CLOSE(UNIT=11) + + +C open statement for output files, removed from line 2449 ff to here + OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000) + GOTO 6001 + 6000 WRITE(*,*)' File schon vorhanden, Gib neue Ausgabedatei an (A8)' + READ(*,'(A8)') fileout + outnam=fileout//outext + rgenam=fileout//rgeext + OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000) + 6001 OPEN(UNIT=22,FILE=rgenam,STATUS='new') + WRITE(21,1000) + 1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - V TrimSP7L 17.Oct.02 TP *') + + CALL TimeStamp(day_start,month_start,year_start,hour_start + & ,min_start,sec_start,seconds_start_total) + WRITE(21,*) + WRITE(21,10050)day_start,month_start,year_start,hour_start + & ,min_start,sec_start + WRITE(*,10050)day_start,month_start,year_start,hour_start + & ,min_start,sec_start +10050 FORMAT(1x,' TrimSp simulation started at: ',A2,'.',A4,1x,A4,1x,A2 + & ,':',A2,':',A2) + +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.) GO TO 2 + E0DE = 100.0D0/(E0*DE) + GO TO 3 + 2 E0DE = 10.0D0/(DABS(E0)*DE) + 3 BW = 180.D0/PI + DAW = BW/DA + DGW = BW/DG + DGIK = BW/DGI + +C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS + PI2=2.D0*PI + ABC=AB*FP + LMAX=7 + JMAX=5 + L=ISRCHEQ(LMAX,DX(1),1,0.D0)-1 + +C Checks wether depth interval is an integer denominator of layer thickness or not +C If not, calculated implantation profile is not correct. + + depth_interval_flag = 1 + DO K=1,L-1 + IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN + depth_interval_flag = 0 + GO TO 44 + ENDIF + ENDDO + 44 CONTINUE + + DO I=1,L + DO J=1,JMAX + IF(EQUAL(CO(I,J),0.D0)) GOTO 156 + ENDDO + J=JMAX+1 + 156 NJ(I)=J-1 + ENDDO + JT(1) = 0 + JT(2) = NJ(1) + JT(3) = NJ(1)+NJ(2) + JT(4) = JT(3)+NJ(3) + JT(5) = JT(4)+NJ(4) + JT(6) = JT(5)+NJ(5) + JT(7) = JT(6)+NJ(6) + LJ = NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+NJ(7) + XX(1)=DX(1) + DO I=2,L + XX(I)=XX(I-1)+DX(I) + ENDDO + DO I=1,L + 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 + ENDDO + DO LL=1,L + ARHO(LL) = RHO(LL)*AN/M2(LL) + LM(LL) = ARHO(LL)**(-1.D0/3.D0) + ASIG(LL) = LM(LL)*ARHO(LL) + PDMAX(LL) = LM(LL)/DSQRT(PI) + K2(LL) = .133743D0*Z2(LL)**(2.D0/3.D0)/DSQRT(M2(LL)) + AM(LL) = CA*ABC*(Z2(LL)**(-1.D0/3.D0)) + FM(LL) = AM(LL)*M2(LL)/(Z1*Z2(LL)*E2*(M1+M2(LL))) + EPS0(LL) = FM(LL)*E0 + ENDDO + DO J = 1,NJ(1) + ZZ(J) = ZT(1,J) + TM(J) = MT(1,J) + DI(J) = ED(1,J) + EP(J) = BE(1,J) + ENDDO + DO J = 1,NJ(2) + ZZ(NJ(1)+J) = ZT(2,J) + TM(NJ(1)+J) = MT(2,J) + DI(NJ(1)+J) = ED(2,J) + EP(NJ(1)+J) = BE(2,J) + ENDDO + DO J = 1,NJ(3) + ZZ(NJ(1)+NJ(2)+J) = ZT(3,J) + TM(NJ(1)+NJ(2)+J) = MT(3,J) + DI(NJ(1)+NJ(2)+J) = ED(3,J) + EP(NJ(1)+NJ(2)+J) = BE(3,J) + ENDDO + DO J = 1,NJ(4) + ZZ(NJ(1)+NJ(2)+NJ(3)+J) = ZT(4,J) + TM(NJ(1)+NJ(2)+NJ(3)+J) = MT(4,J) + DI(NJ(1)+NJ(2)+NJ(3)+J) = ED(4,J) + EP(NJ(1)+NJ(2)+NJ(3)+J) = BE(4,J) + ENDDO + DO J = 1,NJ(5) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ZT(5,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = MT(5,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ED(5,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = BE(5,J) + ENDDO + DO J = 1,NJ(6) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ZT(6,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = MT(6,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ED(6,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = BE(6,J) + ENDDO + DO J = 1,NJ(7) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ZT(7,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = MT(7,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ED(7,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = BE(7,J) + ENDDO + DO I=1,L + 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.3) + 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.1) 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.2) 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.3) 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.3) +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.1) 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.2) 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.3) 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 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 +C MAXWELLIAN VELOCITY DISTRIBUTION +C + CALL VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,NUM) + DO IV=1,NUM + EMX = EMX+E(IV) + ENDDO + DO iv=1,num + ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) + me(ne) = me(ne)+1 + ENDDO + GO TO 56 +C +C MAXWELLIAN ENERGY DISTRIBUTION +C + 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.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) + ENDDO + DO IV=1,NUM + LLL(IV) = JL + ENDDO +C +C PROJECTILE LOOP +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 COLLISION LOOP (INCLUDES WEAK SIMULTANEOUS COLL. FOR KK1.LT.4) +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)) + & +JT(LLL(IV)) + ENDDO + DO IV=1,IH1 + 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.1) 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.2) 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.3) 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.2.OR.KDEE1.EQ.3) + 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 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 + 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)) + ENDDO + DO IV=1,IH1 + 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.0) GO TO 27 +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) + ENDDO + + DO IREC1=1,NREC2 + 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))+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.1) 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.2) 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.3) 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.1) + 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 + 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)))) + ENDDO + DO IREC1=1,NREC2 + 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 +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 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) GO TO 90 + 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), 100), 1) + IPL(IP)=IPL(IP)+1 + I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), 101), 0) + IRP(I1)=IRP(I1)+1 +C +C Berechnung der gestoppten Teilchen im jeweiligen 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) + IPLT(IP)=IPLT(IP)+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 + EMAT(IG,IAGB) = EMAT(IG,IAGB)+ES + EMAT(IG,22) = EMAT(IG,22)+ES + EMAT(NG,IAGB) = EMAT(NG,IAGB)+ES + 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.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 + IF(IVMAX.LT.IVMIN) GO TO 132 + DO IV=IVMAX+1,IH1 + LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) + ENDDO + DO IV=IVMAX+1,IH1 + 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 PRINTOUT +C +C + do ima = 5000,1,-1 + if(me(ima).ne.0) goto 1010 + enddo + ima = 1 + 1010 ima = MIN0(ima+2,5000) + open(20,file='edist') + do ne=1,ima + write(20,1020) ne,me(ne) + enddo + 1020 format(1x,2i6) + close(20) +C +C Berechnung der part. reflec. coeff. nach Thomas et al. +C + E0keV=E0/1.D3 + EsigkeV=Esig/1.D3 + + IF(ZT(1,2).LT.1.0D-3) 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 how many seconds are needed for the simulation ?? +C + CALL TimeStamp(day_stop,month_stop,year_stop, hour_stop,min_stop + & ,sec_stop,seconds_stop_total) + WRITE(21,*) + WRITE(21,10051)day_stop,month_stop,year_stop, hour_stop,min_stop + & ,sec_stop + WRITE(*,10051)day_stop,month_stop,year_stop, hour_stop,min_stop + & ,sec_stop +10051 FORMAT(1x,' TrimSp simulation ended at: ',A2,'.',A4,1x,A4, 1x,A2 + & ,':',A2,':',A2) + WRITE(21,*) + WRITE(21,10052)nh,(seconds_stop_total-seconds_start_total) +10052 FORMAT(1x,' Simulation needed for ',I7,' muons ',I7,' seconds') + + WRITE(21,1402)innam + 1402 FORMAT(//30X,'* INPUT DATA *',5X,A12) + WRITE(21,1404) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC + 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) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT + & ,IPOTR,IRL + 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,L + WRITE(21,1412) I,DX(I),RHO(I),CK(I),(ZT(I,J),J=1,5) ,(MT(I,J),J + & =1,5),(CO(I,J),J=1,5) + 1412 FORMAT(/1X,I1,6H.LAYER,1X,1F8.2,2F7.2,5F6.0,5F7.2,5F6.3) + 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,L + WRITE(21,1418) I,(SBE(I,J),J=1,5),(ED(I,J),J=1,5),(BE(I,J),J=1 + & ,5) + 1418 FORMAT(/1X,I1,6H.LAYER,17X,5F6.2,3X,5F7.2,3X,5F6.2) + 1416 CONTINUE + IF(KDEE1.LT.4) GO TO 1421 + WRITE(21,1419) + 1419 FORMAT(//30X,'CH1',10X,'CH2',10X,'CH3',10X,'CH4',10X,'CH5') + DO 1417 I=1,L + WRITE(21,1415) I,CH1(I,1),CH2(I,1),CH3(I,1),CH4(I,1),CH5(I,1) + IF(NJ(I).LT.2) GO TO 1417 + WRITE(21,1423) CH1(I,2),CH2(I,2),CH3(I,2),CH4(I,2),CH5(I,2) + IF(NJ(I).LT.3) GO TO 1417 + WRITE(21,1423) CH1(I,3),CH2(I,3),CH3(I,3),CH4(I,3),CH5(I,3) + IF(NJ(I).LT.4) GO TO 1417 + WRITE(21,1423) CH1(I,4),CH2(I,4),CH3(I,4),CH4(I,4),CH5(I,4) + IF(NJ(I).LT.5) GO TO 1417 + WRITE(21,1423) CH1(I,5),CH2(I,5),CH3(I,5),CH4(I,5),CH5(I,5) + 1417 CONTINUE + 1415 FORMAT(/1X,I1,6H.LAYER,17X,5F13.6) + 1423 FORMAT(/25X,5F13.6) + 1421 CONTINUE + IF(IPOT.EQ.1) DPOT='KR-C POTENTIAL' + IF(IPOT.EQ.2) DPOT='mod. MOLIERE ' + IF(IPOT.EQ.3) DPOT='ZBL POTENTIAL' + IF(IPOTR.EQ.1) DPOTR='KR-C POTENTIAL' + IF(IPOTR.EQ.2) DPOTR='MOLIERE POTENTIAL' + IF(IPOTR.EQ.3) DPOTR='ZBL POTENTIAL' + WRITE(21,1411) DPOT,DPOTR + 1411 FORMAT(//7X,'INTERACTION POTENTIAL : PROJECTILE-TARGET : ',A18 + & ,' TARGET-TARGET : ',A18) + IF(KDEE1.EQ.1) DKDEE1='LINDHARD-SCHARFF' + IF(KDEE1.EQ.2) DKDEE1='OEN-ROBINSON' + IF(KDEE1.EQ.3) DKDEE1='50% LS 50% OR' + IF(KDEE1.EQ.4) DKDEE1='AZ nach ICRU49' + IF(KDEE1.EQ.5) DKDEE1='ZIEGLER' + IF(KDEE2.EQ.1) DKDEE2='LINDHARD-SCHARFF' + IF(KDEE2.EQ.2) DKDEE2='OEN-ROBINSON' + IF(KDEE2.EQ.3) DKDEE2='50% LS 50% OR' + WRITE(21,1413) DKDEE1,DKDEE2 + 1413 FORMAT(//7X,'INELASTIC LOSS MODEL : PROJECTILE-TARGET : ',A18 + & ,' TARGET-TARGET : ',A18) + IF(E0.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 ALPHAM=-ALPHA + WRITE(21,1407) TI,SHEATH,ALPHAM,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(Esig,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(ALPHASIG,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)outnam + 1428 FORMAT(1H1,//30X,'* OUTPUT DATA *',5X,A12) +#if defined (OS_WIN) + WRITE(22,14280)rgenam +14280 FORMAT(1H1,//30X,'* RANGE DATA *',5X,A12) +#endif + WRITE(21,1430) HLM,HLMT,SU,SUT,XC,RT,SFE,INEL,L,LJ + 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,MAXA,NALL,NPA,NSA,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,L + WRITE(21,473) I,EPS0(I),Z2(I),M2(I),ARHO(I),LM(I),PDMAX(I) + & ,ASIG(I),SB(I),XX(I),NJ(I) + 473 FORMAT(/1X,I1,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,LJ + 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,LJ + 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,LJ + WRITE(21,477) (A(I,J),J=1,LJ) + 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,LJ + WRITE(21,492) (KOR(I,J),J=1,LJ) + 492 FORMAT(/1X,9X,35F8.5) + ENDDO +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) + ENDDO + DO J=1,LJ + ISPAT = ISPAT+ITSP(J) + ESPAT = ESPAT+ETSP(J) + ENDDO + WRITE(21,1500) IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,ISPAT,ESPAT, tryE + & ,negE,epsilon,prcoeff + 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(l.gt.1) then + do i=1,l + nli(i)=IDINT(xx(i)/cw+0.01) + enddo + do i=1,l + do j=nli(i-1)+1,nli(i) + irpl(i)=irpl(i)+irp(j) + enddo + enddo + WRITE(21,15001) +15001 FORMAT(/33x,'FROM BIN WIDTH',2x,'FROM LAYER THICKNESS') + IF(depth_interval_flag.EQ.0) WRITE(21,'(36x,5HWRONG)') + + do i=1,l + WRITE(21,1501) i,irpl(i),number_in_layer(i) + enddo + 1501 FORMAT(/1x,'IMPLANTED PARTICLES (LAYER ',i1,')=',i16,2x,i16) + endif + CSUM=ICSUM + CSUMS=ICSUMS + AVCSUM=CSUM/HN + AVCSMS=CSUMS/HN + AVCDIS=DBLE(ICDI)/HN + CST=DBLE(ICSUM-ICDI) + WRITE(21,1511) AVCSUM,AVCDIS,AVCSMS + 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) + 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 + 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 + WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', + & '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,101 + RIRP(I) = DBLE(IRP(I))/YH + 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,100 + 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,100 + ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) + ENDDO + ENDDO + ENDDO + DO I=0,101 + IIRP=IIRP+IRP(I) + TRIRP=TRIRP+RIRP(I) + ENDDO + DO I=1,100 + 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=100,1,-1 + if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 + enddo + im1=1 + 20 im1=min0(im1+2,100) + DO 11 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,F6.0,2x,I10) + D1=D2 + 11 D2=D2+CW + WRITE(21,604) D2-CW,IRP(101),RIRP(101) + 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,100 + 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=100,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,100) + 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.100) ILD(I)=100 + 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 ',I1,' : ',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=100,1,-1 + if (.not.equal(ionr(im3),0.D0)) go to 31 + enddo + im3=1 + 31 im3=MIN0(im3+2,100) + 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 +C +C BACKSCATTERING +C + IF(IB.EQ.0) GO TO 1512 + WRITE(21,1527) + 1527 FORMAT(1H1,//5X,'BACKSCATTERING OF PROJECTILES'/) + BI=DBLE(IB) + BIL=DBLE(IBL) + RN=BI/HN + RE=EB/(HN*EMV) + EMEANR=RE/RN + EMEAN=EB/BI + AVEB=EMEAN + IF (equal(BI,1.0d0))GO TO 1506 + AVNLB=ENUCLB/BI + VANLB=ENL2B/BI-AVNLB*AVNLB + SIGNLB=DSQRT(VANLB) + DFINLB=SIGNLB/BI + AVILB=EINELB/BI + VAILB=EIL2B/BI-AVILB*AVILB + SIGILB=DSQRT(VAILB) + DFIILB=SIGILB/BI + 1506 WRITE(21,1508) RN,RE,EMEANR,EMEAN + 1508 FORMAT(/5X,'PART.REFL.COEF.=',1PE11.4,' ENERGY REFL.COEF.=' + & ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' + & ,1E11.4) + IF(IB.EQ.0) GO TO 1512 + CALL MOMENT(EB1B,EB2B,EB3B,EB4B,EB5B,EB6B,EB,EB2SUM,EB3SUM,EB4SUM + & ,EB5SUM,EB6SUM,BI) + CALL MOMENT(EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL,EB1SUL,EB2SUL + & ,EB3SUL,EB4SUL,EB5SUL,EB6SUL,BIL) + CALL MOMENTS(FIB0,SEB,THB,FOB,FIB,SIB,SIGMAB,DFIB0,DSEB,DTHB,EB + & ,EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM,BI) + CALL MOMENT(PL1S,PL2S,PL3S,PL4S,PL5S,PL6S,PLSB,PL2SB,PL3SB,PL4SB + & ,PL5SB,PL6SB,BI) + CALL MOMENTS(FIPB0,SEPB,THPB,FOPB,FIPB,SIPB,SIGMPB,DFIPB0,DSEPB + & ,DTHPB, PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB,BI) + WRITE(21,7117) + WRITE(21,7241) FIB0,SEB,THB,FOB,SIGMAB,DFIB0,DSEB,DTHB + 7241 FORMAT(1X,' ENERGY',5X,1P1E12.4,7E14.4) + WRITE(21,7242) FIPB0,SEPB,THPB,FOPB,SIGMPB,DFIPB0,DSEPB,DTHPB + 7242 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + WRITE(21,7237) AVNLB,VANLB,SIGNLB,DFINLB + WRITE(21,7238) AVILB,VAILB,SIGILB,DFIILB + WRITE(21,7118) + WRITE(21,1541) EB1B,EB2B,EB3B,EB4B,EB5B,EB6B + 1541 FORMAT(1X,' ENERGY',5X,1P1E12.4,5E14.4) + WRITE(21,1543) EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL + 1543 FORMAT(1X,' LOGENERGY',5X,1P1E12.4,5E14.4) + WRITE(21,1545) PL1S,PL2S,PL3S,PL4S,PL5S,PL6S + 1545 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + DO 1510 I=1,20 + 1510 AI(I)=0.05D0*DBLE(I) + IF(IB.EQ.0) GO TO 1512 + WRITE(21,1514) + 1514 FORMAT(//5X,'POLAR ANGULAR DISTRIBUTION OF BACKSCATTERED ', + & 'PROJECTILES'//) + DO I=1,20 + RKADB(I)=DBLE(KADB(I))*20.D0/DBLE(IB) + ENDDO + WRITE(21,1518)(AI(I),I=1,20),(KADB(I),I=1,20),(RKADB(I),I=1,20) + 1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) + 1512 CONTINUE +C +C TRANSMISSION +C + IF(IT.EQ.0) GO TO 1524 + WRITE(21,1529) + 1529 FORMAT(///5X,' TRANSMISSION OF PROJECTILES'/) + TIT=DBLE(IT) + TN=TIT/HN + TE=ET/(HN*E0) + TMEANR=TE/TN + EMEANT=TMEANR*E0 + IF (equal(TIT,1.0D0)) GO TO 1520 + AVNLT=ENUCLT/TIT + VANLT=ENL2T/TIT-AVNLT*AVNLT + SIGNLT=DSQRT(VANLT) + DFINLT=SIGNLT/TIT + AVILT=EINELT/TIT + VAILT=EIL2T/TIT-AVILT*AVILT + SIGILT=DSQRT(VAILT) + DFIILT=SIGILT/TIT + 1520 WRITE(21,1522) TN,TE,TMEANR,EMEANT + 1522 FORMAT(//5X,'PART.TRANSM.COEF.=',1PE11.4,' ENERGY TRANSM.COEF.=' + & ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' + & ,1E11.4) + CALL MOMENTS(FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT,ET + & ,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM,TIT) + CALL MOMENTS(FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT,DFIPT0,DSEPT + & ,DTHPT, PLST,PL2ST,PL3ST,PL4ST,PL5ST,PL6ST,TIT) + WRITE(21,7117) + WRITE(21,7241) FIT0,SET,THT,FOT,SIGMAT,DFIT0,DSET,DTHT + WRITE(21,7242) FIPT0,SEPT,THPT,FOPT,SIGMPT,DFIPT0,DSEPT,DTHPT + WRITE(21,7237) AVNLT,VANLT,SIGNLT,DFINLT + WRITE(21,7238) AVILT,VAILT,SIGILT,DFIILT + WRITE(21,1526) + 1526 FORMAT(//5X,'POLAR ANGULAR DISTRIBUTION OF TRANSMITTED ', + & 'PARTICLES'//) + DO I=1,20 + RKADT(I)=DBLE(KADT(I))*20.D0/DBLE(IT) + ENDDO + WRITE(21,1530) (AI(I),I=1,20),(KADT(I),I=1,20),(RKADT(I),I=1,20) + 1530 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) + 1524 CONTINUE +C +C BACKWARD SPUTTERING : YIELDS AND ENERGIES +C + IF(ISPA.EQ.0) GO TO 1700 + 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 (',I1,'.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) GO TO 1700 + 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 +C +C BACKWARD SPUTTERING : ANGULAR DISTRIBUTIONS +C + WRITE(21,1601) + 1601 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF ALL BACKWARD ', + & 'SPUTTERED PARTICLES'//) + DO 1603 I=1,20 + 1603 RKADS(I)=KADS(I)*20.D0/ISPA + 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 1608 I=1,20 + 1608 RKADSL(I,1)=KADSL(I,1)*20.D0/ISPAL(1) + 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 1620 I=1,20 + 1620 RKADSJ(I,J)=KADSJ(I,J)*20.D0/IBSP(J) + 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 1612 I=1,20 + 1612 RKADSL(I,2)=KADSL(I,2)*20.D0/ISPAL(2) + 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 1628 I=1,20 + 1628 RKADSJ(I,J)=KADSJ(I,J)*20.D0/IBSP(J) + WRITE(21,1518) (AI(I),I=1,20),(KADSJ(I,J),I=1,20) + & ,(RKADSJ(I,J),I=1,20) + 1624 CONTINUE + 1622 CONTINUE +C +C TRANSMISSION SPUTTERING : YIELDS AND ENERGIES +C + 1700 IF(ISPAT.EQ.0) GO TO 1800 + 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,LJ + ISPALT(3) = ISPALT(3)+ITSP(J) + ESPALT(3) = ESPALT(3)+ETSP(J) + ENDDO + WRITE(21,1712) ISPAT,ESPAT + 1712 FORMAT(///,8X,'ALL SPUTTERED PARTICLES = ',I7,3X + & ,'TOTAL SPUTTERED ENERGY = ',E10.4,3H EV//) + DO J=1,L + WRITE(21,1713) J,ISPALT(J),ESPALT(J) + 1713 FORMAT(8X,'SPUTTERED PARTICLES (LAYER ',I1,') = ',I7,3X + & ,'SPUTTERED ENERGY = ',E10.4,3H EV) + ENDDO + DO J=1,LJ + 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)/ESPAT + REIST(J)=ESPIST(J)/ESPAT + REOPT(J)=ESPOPT(J)/ESPAT + REOST(J)=ESPOST(J)/ESPAT + ENDDO + 1715 CONTINUE + DO 1717 J=1,LJ + 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,LJ + SPYT(J)=DBLE(ITSP(J))/DBLE(NH) + SPET(J)=ETSP(J)/(NH*E0) + ENDDO + DO 1737 J=1,LJ + IF (equal(SPYT(J),0.0D0))GO TO 1737 + REYT(J)=SPET(J)/SPYT(J) + EMSPT(J)=REYT(J)*E0 + 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(L.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(L.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,LJ + WRITE(21,1564) J-JT(3),ITSP(J),J-JT(3),ETSP(J) + ENDDO + DO J=JT(3)+1,LJ + 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,LJ + 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 +C +C TRANSMISSION SPUTTERING : ANGULAR DISTRIBUTIONS +C + 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) (AI(I),I=1,20),(KADST(I),I=1,20),(RKADST(I),I=1,20) + IF(L.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 J=NJ(1)+1,JT(3) + KDSTL(I,2)=KDSTL(I,2)+KDSTJ(I,J) + 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) (AI(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) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) + & ,(RKDSTJ(I,J),I=1,20) + 1780 CONTINUE + 1778 IF(L.EQ.1) GO TO 1800 + WRITE(21,1786) + 1786 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ', + & 'PARTICLES ; LAYER 2'//) + WRITE(21,1518) (AI(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) (AI(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,LJ-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,LJ-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) (AI(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) (AI(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) (AI(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,LJ-NJ(1) + WRITE(21,1783) J-NJ(2) + 1783 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED ', + & 'PARTICLES ; LAYER 3 ; SPECIES ',I1//) + WRITE(21,1518) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) + & ,(RKDSTJ(I,J),I=1,20) + ENDDO + 1800 CONTINUE +c +c The file for33 is created here +c + DO i=1,100 + READ(33,'(A246)',ERR=7800,END=7800)COLUMN(i) + ENDDO + 7800 COLCOUNT=i-1 + CLOSE(33,STATUS='DELETE') + WRITE(33,7802) + 7802 FORMAT(6x,'Energy',4x,'SigmaE',5x,'Alpha',2x,'SigAlpha',4x,'ntot', + & 5x,'imp',2x,'backsc',3x,'trans',3x,'tried',4x,'negE',3x, + & 'impL1',3x,'impL2',3x,'impL3',3x,'impL4',3x,'impL5',3x + & ,'impL6', 3x,'impL7',3x, 'range',6x,'straggeling',2x, 'Eback' + & ,7x,'sigEback',4x,'Etrans',6x,'SigEtrans',3x, 'red. E',6x + & ,'PRC') + DO i=2,COLCOUNT + WRITE(33,'(A246)')COLUMN(i) + ENDDO + IF(l.EQ.1) THEN + number_in_layer(1)=IIM + DO k=2,7 + number_in_layer(k)=0 + ENDDO + ELSEIF(l.EQ.2) THEN + DO k=3,7 + number_in_layer(k)=0 + ENDDO + ELSEIF(l.EQ.3) THEN + DO k=4,7 + number_in_layer(k)=0 + ENDDO + ELSEIF(l.EQ.4) THEN + DO k=5,7 + number_in_layer(k)=0 + ENDDO + ELSEIF(l.EQ.5) THEN + DO k=6,7 + number_in_layer(k)=0 + ENDDO + ELSEIF(l.EQ.6) THEN + number_in_layer(7)=0 + ENDIF + WRITE(33,7801)E0keV,EsigkeV,ALPHA,ALPHASIG,NH,IIM,IB,IT,tryE,negE + & ,(number_in_layer(k),k=1,7),FIX0,SIGMAX,FIB0,SIGMAB,FIT0 + & ,SIGMAT,epsilon,prcoeff + 7801 FORMAT(F12.2,3(1x,F9.2),1x,13(I7,1x),6(E12.4),2(E12.4)) + CLOSE(33) +c +c End of file for33 +C +C TOP AND FRONT LINES FOR MATRICES +C + JE=DE + JA=DA + JG=DG + DO J=2,NG1 + MAGB(J,1) = (J-1)*JG + MAGT(J,1) = (J-1)*JG + EMA(J,1)=DBLE(J-1)*DG + EMAT(J,1)=DBLE(J-1)*DG + 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,101 + 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,NG1 + MAGS(J,1,K) = (J-1)*JG + MAGST(J,1,K) = (J-1)*JG + MAGSA(J,1,K) = (J-1)*JG + ENDDO + DO J=2,NA1 + MAGSA(1,J,K) = (J-1)*JA + 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,101 + 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 IG2=1,NGIK,1 + DO J=2,21 + MEAGS(1,IG2,J,K) = J-1 + ENDDO + DO J=2,101 + MEAGS(J,IG2,1,K) = J-1 + ENDDO + ENDDO + ENDDO + DO IG2=1,NGIK,1 + DO J=2,21 + MEAGB(1,IG2,J) = J-1 + MEAGT(1,IG2,J) = J-1 + ENDDO + DO J=2,101 + MEAGB(J,IG2,1) = J-1 + MEAGT(J,IG2,1) = J-1 + ENDDO + ENDDO + DO I=2,74 + ELOG(I)=10.D0**(I/12.D0)*10.D0**(-7.D0/6.D0) + ENDDO + TEMP=(10.D0**(1.D0/12.D0)-1.D0)*10.D0**(-7.D0/6.D0) + TEMPNH=TEMP*DBLE(NH) +C +C BACKWARD SPUTTERING : MATRICES , ENERGY - ANGLE CORRELATIONS +C + IF(ISPA.LT.10000) GO TO 1900 + DO 1850 J=1,JT(3) + EASL(2,J)=DBLE(MEASL(2,21,J))/(DBLE(NH)*0.1) + DO 1850 IESLOG=3,74 + 1850 EASL(IESLOG,J)=DBLE(MEASL(IESLOG,21,J))/(TEMPNH*10.D0 + & **((IESLOG-1)/12.D0)) + 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 = 101,1,-1 + if(meas(ima,22,j).ne.0) goto 1883 + enddo + ima = 1 + 1883 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meags(ima,ig2,22,j).ne.0) goto 1885 + enddo + ima = 1 + 1885 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meas(ima,22,j).ne.0) goto 1895 + enddo + ima = 1 + 1895 ima = min(ima+2,101) + 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 +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 1950 J=1,JTJ + EASTL(2,J)=DBLE(MEASTL(2,21,J))/(DBLE(NH)*0.1D0) + DO 1950 IESLOG=3,74 + 1950 EASTL(IESLOG,J)=DBLE(MEASTL(IESLOG,21,J))/(TEMPNH*10.D0 + & **((IESLOG-1)/12.D0)) + 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 = 101,1,-1 + if(meast(ima,22,j).ne.0) goto 1983 + enddo + ima = 1 + 1983 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meast(ima,22,j).ne.0) goto 1995 + enddo + ima = 1 + 1995 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meast(ima,22,j).ne.0) goto 1977 + enddo + ima = 1 + 1977 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meagb(ima,ig2,22).ne.0) goto 2015 + enddo + ima = 1 + 2015 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meab(ima,22).ne.0) goto 2017 + enddo + ima = 1 + 2017 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meagt(ima,ig2,22).ne.0) goto 2115 + enddo + ima = 1 + 2115 ima = min(ima+2,101) + 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 = 101,1,-1 + if(meat(ima,22).ne.0) goto 2117 + enddo + ima = 1 + 2117 ima = min(ima+2,101) + 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 STOP + END + + + + 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 + + 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 + + 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 + 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 + IMPLICIT NONE + INTEGER INIV1,INIV3 + REAL*8 FG(128),FFG(64) + 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 + 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, 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)) + 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) + 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) + REAL*8 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 + + 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 + + INTEGER FUNCTION ISRCHEQ(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).EQ.TARGT) GO TO 200 + J=J+INC + ENDDO + 200 ISRCHEQ=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/ + 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/ + 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) + 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 + + + INTEGER FUNCTION OldNew(filename) +C This funnction returns 0 for old input format and 1 for new input +C format + CHARACTER filename*12 + REAL dummy(21) + INTEGER idummy +C Assume old format + OldNew = 0 + +C Read first three lines, the third line differs between the two +C format + OPEN(UNIT=11,file=filename,STATUS='unknown') + READ(11,*) + READ(11,*) + READ(11,*,ERR=10) dummy(1),dummy(2),dummy(3),dummy(4),dummy(5) + & ,dummy(6),dummy(7),dummy(8),dummy(9),dummy(10),dummy(11) + & ,dummy(12),dummy(13),dummy(14),dummy(15),dummy(16),dummy(17) + & ,dummy(18),dummy(19),dummy(20),dummy(21) + GOTO 20 +C If you got here it means that the file in in the old format + 10 OldNew = 1 + 20 CLOSE(11) + if (OldNew.eq.1) then + write(*,*) "The input file is in the new format" + else + write(*,*) "The input file is in the old format" + endif + + RETURN + END + + 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 195 I= 1, 24 + NEXT(I) = I-1 + 195 TWOM24 = TWOM24 * 0.5 + 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 +