C Version TrimSp7L ----> 7Layer C C erstellt 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 CDIR$ NOLIST C !DEC$REAL:8 C C IMPLICIT INTEGER (i-j) C IMPLICIT REAL*8 (a-h,k-z) C 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 Date_time(8) INTEGER depth_interval_flag INTEGER*4 days_start_total,days_stop_total 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:1001),IPL(1000),IPLB(1000),IPLT(1000) INTEGER*4 ICD(1000,35),ICDT(1000),ICDJT(35),ICDIRJ(35,35) # ,ICDR(1000,35),ICDTR(1000),ICDJTR(35) # ,ICDIRI(1000,35,35),ICDIRN(1000,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,1002) INTEGER*4 MEAT(102,22),MAGT(62,22),MEAGT(102,36,22), # MEATL(75,21),MEPT(102,1002) 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) CC # ,MEAGST(102,36,22,10) von Eckstein herauskommentiert # ,MEASTL(75,21,30) CC REAL*8 MEAGSL(75,36,21),EAGSL(75) von Eckstein herauskommentiert 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 C REAL Variablen 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:1001) # ,CASMOT(1000),PHON(1000),DENT(1000),ION(1000),DMGN(1000) # ,CASMOTR(1000),PHONR(1000),DENTR(1000),IONR(1000) # ,DMGNR(1000),ELGD(1000),ELGDR(1000) REAL*8 ELE(1000,35),ELI(1000,35),ELP(1000,35),ELD(1000,35) # ,ELET(35),ELIT(35),ELPT(35),ELDT(35) # ,ELER(1000,35),ELIR(1000,35),ELPR(1000,35),ELDR(1000,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) CC REAL*8 SL(64,5),SM(64,5),SH(64,5),CH(92,12),CHE(92,9) 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) CC REAL*8 ESVDL(2000) 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 C CHARACTER Variablen 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 Real_Clock(3)*12 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 C COMMON /A/ M1,VELC,ZARG COMMON /B/ TI,SHEATH,CALFA C DATA PI/3.14159265358979D0/, ICW/100/, E2/14.399651D0/ DATA AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/ DATA inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/ DATA errext/'.err'/ 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/102204*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/1000*0.D0/,DENT/1000*0.D0/ DATA DMGN/1000*0.D0/,ION/1000*0.D0/,PHON/1000*0.D0/ DATA PHONR/1000*0.D0/ DATA ELGD/1000*0.D0/,ELGDR/1000*0.D0/ DATA ICDT/1000*0/,ICDTR/1000*0/ DATA ICDR/35000*0/,ICDIRN/35000*0/,IONR/1000*0.D0/ DATA DENTR/1000*0.D0/,DMGNR/1000*0.D0/ DATA IPL/1000*0/,IPLB/1000*0/,IPLT/1000*0/ c DATA IRP/102*0/ gibt witzigweise einen Fehler, aber warum???? DATA IRPL/7*0/ DATA ICDJT/35*0/,ICDJTR/35*0/,ICDITR/35*0/ DATA ICD/35000*0/,ELP/35000*0.D0/,ELD/35000*0.D0/ DATA ELE/35000*0.D0/,ELI/35000*0.D0/ DATA ICDIRI/1225000*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/ C CC EXTERNAL CVMGT,ILLZ,SCOPY,ISRCHEQ,ISRCHFGE,ISRCHFGT C innam=filein//inext outnam=fileout//outext rgenam=fileout//rgeext errnam=fileout//errext C OPEN(UNIT=99,file=errnam,STATUS='new') OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=13591) C OPEN(UNIT=31,NAME='energ.dat',STATUS='new') C READ(11,100) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC READ(11,101) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 # ,IPOT,IPOTR,IRL READ(11,102) 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) DO 135 I=1,7 READ(11,103) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) READ(11,103) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) READ(11,103) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) READ(11,103) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) READ(11,103) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) READ(11,103) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) READ(11,107) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) READ(11,107) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) READ(11,107) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) READ(11,107) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) READ(11,107) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) 135 CONTINUE 13591 CLOSE(UNIT=21) WRITE(*,100) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC WRITE(*,101) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 # ,IPOT,IPOTR,IRL WRITE(*,102) 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) DO 1359 I=1,7 WRITE(*,'(1x,I2,A7)')i,'. Layer' WRITE(*,103) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) WRITE(*,103) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) WRITE(*,103) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) c WRITE(*,103) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) c WRITE(*,103) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) c WRITE(*,103) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) c WRITE(*,107) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) c WRITE(*,107) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) c WRITE(*,107) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) c WRITE(*,107) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) c WRITE(*,107) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) 1359 CONTINUE C 100 FORMAT(2F7.2,1F12.2,7F9.2) 101 FORMAT(I9,3F8.0,1F7.2,1F7.0,2F7.2,6I4,I3) 102 FORMAT(7F9.2,14F7.2) 103 FORMAT(5F9.4) 107 FORMAT(5F12.6) C C open statement for output files, removed from line 2449 ff to here C 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 26.06.00 **') C C 1st CALL DATE_AND_TIME CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3), 1 Date_Time) C IF(Date_Time(2).EQ.1) THEN month_start='Jan.' days_start_total=Date_Time(3) ELSEIF(Date_Time(2).EQ.2) THEN month_start='Feb.' days_start_total=31+Date_Time(3) ELSEIF(Date_Time(2).EQ.3) THEN month_start='Mar.' days_start_total=31+28+Date_Time(3) ELSEIF(Date_Time(2).EQ.4) THEN month_start='Apr.' days_start_total=31+28+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.5) THEN month_start='May ' days_start_total=31+28+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.6) THEN month_start='Jun.' days_start_total=31+28+31+30+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.7) THEN month_start='Jul.' days_start_total=31+28+31+30+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.8) THEN month_start='Aug.' days_start_total=31+28+31+30+31+30+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.9) THEN month_start='Sep.' days_start_total=31+28+31+30+31+30+31+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.10) THEN month_start='Oct.' days_start_total=31+28+31+30+31+30+31+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.11) THEN month_start='Nov.' days_start_total=31+28+31+30+31+30+31+31+30+31+Date_Time(3) ELSE month_start='Dec.' days_start_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3) ENDIF C in seconds from beginning of year seconds_start_total=Date_Time(7)+(Date_Time(6)*60)+ 1 (Date_Time(5)*3600)+(days_start_total-1)*86400 C READ(Real_Clock(1)(1:4),'(A4)')year_start READ(Real_Clock(1)(7:8),'(A2)')day_start READ(Real_Clock(2)(1:2),'(A2)')hour_start READ(Real_Clock(2)(3:4),'(A2)')min_start READ(Real_Clock(2)(5:6),'(A2)')sec_start C WRITE(21,*) WRITE(21,10050)day_start,month_start,year_start, 1 hour_start,min_start,sec_start 10050 FORMAT(1x,' TrimSp simulation started at: ',A2,'.',A4,1x,A4, 1 1x,A2,':',A2,':',A2) C SET INTERVAL CONSTANTS FOR OUTPUT C 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 C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS C PI2=2.D0*PI ABC=AB*FP LMAX=7 JMAX=5 L=ISRCHEQ(LMAX,DX(1),1,0.D0)-1 C C Checks wether depth interval is an integer denominator of layer thickness or not C If not, calculated implantation profile is not correct. C depth_interval_flag = 1 LOOP_Check_layer_thick : DO K=1,L-1 IF(.NOT.EQUAL(DX(K)/CW-DFLOAT(JFIX(DX(K)/CW)),0.D0)) THEN depth_interval_flag = 0 EXIT LOOP_Check_layer_thick ENDIF ENDDO LOOP_Check_layer_thick C DO 165 I=1,L DO 155 J=1,JMAX IF(EQUAL(CO(I,J),0.D0)) GOTO 156 C IF(CO(I,J).D0EQ.D00.D000) GO TO 156 155 CONTINUE J=JMAX+1 156 NJ(I)=J-1 165 CONTINUE 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 170 I=2,L 170 XX(I)=XX(I-1)+DX(I) DO 180 I=1,L DO 180 J=1,NJ(I) Z2(I)=Z2(I)+CO(I,J)*ZT(I,J) M2(I)=M2(I)+CO(I,J)*MT(I,J) 180 CONTINUE DO 185 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 185 CONTINUE DO 186 J = 1,NJ(1) ZZ(J) = ZT(1,J) TM(J) = MT(1,J) DI(J) = ED(1,J) 186 EP(J) = BE(1,J) DO 187 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) 187 EP(NJ(1)+J) = BE(2,J) DO 188 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) 188 EP(NJ(1)+NJ(2)+J) = BE(3,J) DO 1880 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) 1880 EP(NJ(1)+NJ(2)+NJ(3)+J) = BE(4,J) DO 1881 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) 1881 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = BE(5,J) DO 1882 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) 1882 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = BE(6,J) DO 18803 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) 18803 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = BE(7,J) DO 189 I=1,L COM(1,I) = CO(I,1) DO 189 J=1,NJ(I)-1 COM(J+1,I) = COM(J,I)+CO(I,J+1) 189 CONTINUE DO 190 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)) 190 CONTINUE IF(IPOT.EQ.1) THEN C KR-C POTENTIAL (IPOT=1) DO 194 J=1,LJ KOR1(J) = 0.0389205D0*KL1(J)/(PI*A1(J)*A1(J)) 194 CONTINUE ELSEIF (IPOT.EQ.2) THEN C MOLIERE POTENTIAL (IPOT=2) DO 195 J=1,LJ KOR1(J) = 0.045D0*KL1(J)/(PI*A1(J)*A1(J)) 195 CONTINUE ELSEIF (IPOT.EQ.3) THEN C ZBL POTENTIAL DO 196 J=1,LJ KOR1(J) = 0.0203253D0*KL1(J)/(PI*A1(J)*A1(J)) 196 CONTINUE ENDIF DO 191 I = 1,LJ DO 191 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))) 191 CONTINUE IF (IPOTR.EQ.1) THEN C KR-C POTENTIAL (IPOTR=1) DO 197 I = 1,LJ DO 197 J = 1,LJ KOR(I,J) = 0.0389205D0*KL(I,J)/(PI*A(I,J)*A(I,J)) 197 CONTINUE ELSEIF (IPOTR.EQ.2) THEN C MOLIERE POTENTIAL (IPOTR=2) DO 198 I = 1,LJ DO 198 J = 1,LJ KOR(I,J) = 0.045D0*KL(I,J)/(PI*A(I,J)*A(I,J)) 198 CONTINUE ELSEIF (IPOTR.EQ.3) THEN C ZBL POTENTIAL (IPOTR=3) DO 184 I = 1,LJ DO 184 J = 1,LJ KOR(I,J) = 0.0203253D0*KL(I,J)/(PI*A(I,J)*A(I,J)) 184 CONTINUE ENDIF DO 192 LL=1,L DO 192 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) 192 CONTINUE DO 193 I=1,LJ DO 193 LL = 1,L DO 193 J=1,NJ(LL) C KLM(LL,I) = KLM(LL,I)+CK(LL)*CO(LL,J)*KL(I,J+JT(LL)) KLM(LL,I) = KLM(LL,I)+CO(LL,J)*KL(I,J+JT(LL)) 193 CONTINUE C C ALPHA = CVMGT( .001, ALPHA, ALPHA.EQ.0. ) C ALPHA = CVMGT( 179.999, ALPHA, ALPHA.EQ.180.) 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 C SET CONSTANT DISTANCES C 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 C 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) CC ANFANG = RANSET(IY) CC ANFANG = SRAND48(IY) ISEED = IY ISEED2 = IY2 ISEED3 = IY3 C WRITE(*,*)ISEED C 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 49 IV=1,NUM EMX = EMX+E(IV) 49 CONTINUE DO iv=1,num ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) me(ne) = me(ne)+1 ENDDO c GO TO 56 C C MAXWELLIAN ENERGY DISTRIBUTION C ELSE C CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM) DO 48 IV=1,NUM EMX = EMX+E(IV) 48 CONTINUE GO TO 56 ENDIF C 47 CONTINUE C IF(EQUAL(Esig,0.D0)) THEN C FIXED PROJECTILE ENERGY DO IV=1,NUM E(IV) = E0 C WRITE(*,*)' Da Esig=0 ist E=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 C WRITE(*,*)E(IV),Epar,E0 ENDDO ENDIF C C die nachfolgende Zeile wurden von Linie 633 hier hin verschoben C SFE = DMIN1(SB(1),SB(L)) C IF ( ALPHA.GE.0.D0 ) THEN C IF(EQUAL(ALPHASIG,0.D0))THEN C C FIXED PROJECTILE ANGLE C C nachfolgende drei Zeilen waren vorher bei LINIE 633 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) C WRITE(88,*)ALPHA,ALPHASIG,CALFA,SALFA 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) C WRITE(88,'(5(F12.5))')ALPHA,ALPHASIG, C + ALPHApar,CALFA,SALFA ENDDO ENDIF C ELSEIF (EQUAL(ALPHA,-2.D0)) THEN C ELSEIF ( ALPHA.EQ.-2. ) THEN C C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL) C CDIR$ IVDEP DO 53 IV=1,NUM CC RPHI = PI2*RANF() CC RPHI = PI2*DRAND48() RPHI = PI2*DBLE(RAN(ISEED)) CC RTHETA = RANF() CC RTHETA = DRAND48() RTHETA = DBLE(RAN(ISEED)) COSX(IV) = DSQRT(RTHETA) SINE(IV) = DSQRT(1.D0-RTHETA) COSY(IV) = SINE(IV)*DCOS(RPHI) COSZ(IV) = SINE(IV)*DSIN(RPHI) 53 CONTINUE C ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN C ELSEIF ( ALPHA.EQ.-1. AND. X0.GT.0. ) THEN C C RANDOM DISTRIBUTION C CDIR$ IVDEP DO 50 IV=1,NUM CC RPHI = PI2*RANF() CC RPHI = PI2*DRAND48() RPHI = PI2*DBLE(RAN(ISEED)) CC RTHETA = RANF() CC RTHETA = DRAND48() RTHETA = DBLE(RAN(ISEED)) 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) 50 CONTINUE C C ELSEIF ( ALPHA.EQ.-1. AND. X0.LE.0. ) THEN ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.LE.0.D0) THEN C CDIR$ IVDEP DO 55 IV=1,NUM CC RPHI = PI2*RANF() CC RPHI = PI2*DRAND48() RPHI = PI2*DBLE(RAN(ISEED)) CC RTHETA = RANF() CC RTHETA = DRAND48() RTHETA = DBLE(RAN(ISEED)) 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) 55 CONTINUE C ENDIF C 56 IF ( X0.GT.0.D0 ) GO TO 59 C C EXTERNAL START C DO 57 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 57 CONTINUE C C LOCUS OF FIRST COLLISION C 59 JL = ISRCHFGT(L,XX(1),1,X0) C WRITE(*,*)X0 DO 58 IV=1,NUM CC RA = CVMGT(RANF(),1.0,X0.LE.0.0) CC RA = CVMGT(DRAND48(),1.0,X0.LE.0.0) RA = CVMGT(DBLE(RAN(ISEED)),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) 58 CONTINUE C DO 199 IV=1,NUM LLL(IV) = JL 199 CONTINUE C C PROJECTILE LOOP C 1 CONTINUE C NPROJ=NPROJ+1 C DO 63 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 63 CONTINUE 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 298 IV=1,IH1 JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1 CC # ,RANF())+JT(LLL(IV)) CC # ,DRAND48())+JT(LLL(IV)) # ,DBLE(RAN(ISEED)))+JT(LLL(IV)) 298 CONTINUE DO 67 IV=1,IH1 EPS(IV)=E(IV)*F1(JJJ(IV)) 67 CONTINUE C CDIR$ IVDEP DO 64 IV=1,IH1 C C RANDOM AZIMUTHAL ANGLE AND IMPACT PARAMETER C CC PHIP=PI2*RANF() CC PHIP=PI2*DRAND48() PHIP=PI2*DBLE(RAN(ISEED)) CPHI(IV)=DCOS(PHIP) SPHI(IV)=DSIN(PHIP) CC P(IV)=PDMAX(LLL(IV))*DSQRT(RANF()+KK) CC P(IV)=PDMAX(LLL(IV))*DSQRT(DRAND48()+KK) P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(RAN(ISEED))+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) C IF(A1(JJJ(IV)).EQ.0.) WRITE(99,'(A50)')' A1 vor Label 64 ' B(IV)=P(IV)/A1(JJJ(IV)) 64 CONTINUE CALL SCOPY(IH1,B,1,R,1) C WRITE(99,*)IH1,B(IV),R(IV) C CALL MAGICKRC(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 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 C CALL MAGICKRC(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 104 DO 105 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) C IF(R(IV).EQ.0.D0)WRITE(99,'(1x,F15.7,1x,A50)')R(IV),' Label 104 ' 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)+ 1 .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 105 CONTINUE 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 108 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) 108 S2(IV)=1.D0-(1.D0*C2(IV)) GO TO 4103 C 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 4105 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)) C IF(R(IV).EQ.0.D0)WRITE(99,'(A50)')' R nach Label 4104 ' 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 4105 CONTINUE 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 4108 IV=1,IH1 C IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 4106 ' 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) 4108 S2(IV)=1.D0-(1.D0*C2(IV)) GO TO 4103 C 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 5105 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)) C IF(R(IV).EQ.0.D0)WRITE(99,'(A50)')' R nach Label 5104 ' RR1=1.D0/R(IV) V(IV)=(.02817D0*EX1(IV)+.28022D0*EX2(IV)+.50986D0*EX3(IV)+ 1 .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)+ 1 .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 5105 CONTINUE 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 5108 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) 5108 S2(IV)=1.D0-(1.D0*C2(IV)) 4103 CONTINUE C C END OF MAGIC C DO 65 IV=1,IH1 DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV) C TAU(IV)=CVMGT(P(IV)*DSQRT(S2(IV)/C2(IV)),0.,KK.EQ.4) IF(C2(IV).LT.1.D-10) THEN c WRITE(*,*)C2(IV),S2(IV) 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) c WRITE(*,*)C2(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 65 CONTINUE 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 14 IV=1,IH1 ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV)) TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV)) 14 CONTINUE GO TO(15,16,17,18,19),KDEE1 15 DO 26 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) 26 CONTINUE GO TO 40 16 DO 21 IV=1,IH1 DEE(IV)=DEES(IV) 21 CONTINUE GO TO 40 17 DO 22 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) 22 CONTINUE GO TO 40 18 DO 23 IV=1,IH1 SM(IV)=0.D0 EM(IV)=E(IV)*0.001D0/M1 23 CONTINUE DO 66 IV=1,IH1 DO 66 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) 66 CONTINUE DO 73 IV=1,IH1 DO 73 J=1,NJ(LLL(IV)) SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) 73 CONTINUE DO 78 IV=1,IH1 DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10.D0) 78 CONTINUE DO 69 IV=1,IH1 DEE(IV)=10.D0*ASIGT(IV)* # CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT) 69 CONTINUE GO TO 40 19 FHE=CVMGT(1.3333D0,1.D0,M1.LT.4.00D0) DO 25 IV=1,IH1 SM(IV)=0.D0 EM(IV)=E(IV)*0.001D0*FHE 25 CONTINUE DO 74 IV=1,IH1 DO 74 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))) 74 CONTINUE DO 92 IV=1,IH1 DO 92 J=1,NJ(LLL(IV)) SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) 92 CONTINUE DO 79 IV=1,IH1 DEE(IV)=10.D0*ASIGT(IV)* # CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT) 79 CONTINUE 40 CONTINUE C DO 44 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)) 44 CONTINUE C C INCREMENT OF DAMAGE, CASCADE AND PHONON ENERGY C DO 70 IV=1,IH1 C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 I=MAX0(MIN0(JFIX(X1(IV)/CW+1.D0),1000),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 80 IV=1,IH1 ICDI=ICDI+JFIX(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV)))) ICSUMS=ICSUMS+JFIX(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1))) ICSUM=ICSUM+JFIX(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0)) 80 CONTINUE DO 72 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) 72 CONTINUE IF(KK0.EQ.0) GO TO 89 DO 71 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 71 CONTINUE 89 CONTINUE C C IF IRL=0 NO RECOILS ARE FOLLOWED 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 cc IF(DEN(IV).LE.SFE) GO TO 6 IF(DEN(IV).LE.ERC) GO TO 6 IF(X1(IV).GT.RD.AND.X1(IV).LT.RT) GO TO 6 C C CALL NEWREC(NREC1,DEN(IV),X(IV),Y(IV),Z(IV), C 1 CX(IV),CY(IV),CZ(IV),SX(IV), C 2 CT(IV),ST(IV),PHI(IV),P(IV), C 3 ER(1,1),XR(1,1),YR(1,1),ZR(1,1),PHIR(1,1),PSIR(1,1), C 4 CSXR(1,1),CSYR(1,1),CSZR(1,1),SNXR(1,1),L(1,1) 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 CDIR$ IVDEP DO 91 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) 91 CONTINUE 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 1 INCREASED') 8886 CONTINUE C C PROCESS THE PARTICLES IN LIST 2 C C FIND LAYER C DO 68 IREC1=1,NREC2 LRR(IREC1,2)=MIN0(ISRCHFGT(L,XX(1),1,XR(IREC1,2)),L) 68 CONTINUE C C MOVE PARTICLES C DO 60 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) 60 CONTINUE C DO 81 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 81 CONTINUE C KK2=KK0R DO 235 KKR=KK2,0,-1 C C CHOICE OF COLLISION PARTNERS C DO 303 IREC1=1,NREC2 JJR(IREC1,1) = ISRCHFGE(NJ(LRR(IREC1,2)),COM(1,LRR(IREC1,2)) CC 1 ,1,RANF())+JT(LRR(IREC1,2)) CC 1 ,1,DRAND48())+JT(LRR(IREC1,2)) 1 ,1,DBLE(RAN(ISEED)))+JT(LRR(IREC1,2)) 303 CONTINUE C CDIR$ IVDEP DO 236 IREC1=1,NREC2 CC PHIPR=PI2*RANF() CC PHIPR=PI2*DRAND48() PHIPR=PI2*DBLE(RAN(ISEED)) CPHIR(IREC1,2)=DCOS(PHIPR) SPHIR(IREC1,2)=DSIN(PHIPR) CC PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(RANF()+KKR) CC PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DRAND48()+KKR) PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(RAN(ISEED))+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 1 .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) 236 CONTINUE C 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 206 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))/ 1 EPSR(IV)-1.D0 Q=FR/FR1 RR(IV)=RR(IV)-Q TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 206 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 207 IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 207 GO TO 205 207 DO 208 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) 208 S2R(IV)=1.D0-C2R(IV) GO TO 4203 C 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 4206 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))/ 1 EPSR(IV)-1.D0 Q=FR/FR1 RR(IV)=RR(IV)-Q TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 4206 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 4207 IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 4207 GO TO 4205 4207 DO 4208 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) 4208 S2R(IV)=1.D0-C2R(IV) GO TO 4203 C 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)+ 1 .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)+ 1 .4804359794D0*EX3R(IV)+.581563650D0*EX4R(IV))*RRR1 FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ 1 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 5208 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) 5208 S2R(IV)=1.D0-C2R(IV) 4203 CONTINUE C DO 237 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, 1 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 237 CONTINUE C CALL DIRCOS(CSXR(1,2),CSYR(1,2),CSZR(1,2),SNXR(1,2), 1 CPSIR(1,2),SPSIR(1,2),CPHIR(1,2),SPHIR(1,2),NREC2) C 235 CONTINUE C C CREATE SECONDARY KNOCK-ON ATOMS C DO 246 IREC1=1,NREC2 cc IF(T(IREC1).LE.SFE) GO TO 246 IF(T(IREC1).LE.ERC) GO TO 246 IF(X2(IREC1).GT.RD.AND.X2(IREC1).LT.RT) GO TO 246 C C CALL NEWREC(NREC1,T(IREC1),XR(IREC1,2),YR(IREC1,2),ZR(IREC1,2), C 1 CXR(IREC1),CYR(IREC1),CZR(IREC1),SXR(IREC1), C 2 CTR(IREC1),STR(IREC1),PHIR(IREC1,2),PR(IREC1), C 3 ER(I,1),XR(1,1),YR(1,1),ZR(1,1),PHIR(1,1),PSIR(1,1) C 4 ,CSXR(1,1),CSYR(1,1),CSZR(1,1),SNXR(1,1),L(1,1) 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)- 1 CPHIR(IREC1,2)*CYR(IREC1)*CXR(IREC1))/SXR(IREC1) ZR(NREC1,1)=ZR(IREC1,2)+PR(IREC1)*(SPHIR(IREC1,2)*CYR(IREC1)+ 1 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 238 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)) 238 CONTINUE GO TO(115,116,117),KDEE2 115 DO 241 IREC1=1,NREC2 DEER(IREC1)=CVMGT(0.,KLM(LRR(IREC1,2),JJR(IREC1,2))*ASIGTR(IREC1)* # DSQRT(ER(IREC1,2)),XR(IREC1,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) 241 CONTINUE GO TO 242 116 DO 243 IREC1=1,NREC2 DEER(IREC1)=DEERS(IREC1) 243 CONTINUE GO TO 242 117 DO 244 IREC1=1,NREC2 DEER(IREC1)=CVMGT(DEERS(IREC1),.5*(KLM(LRR(IREC1,2),JJR(IREC1,2))* 1 ASIGTR(IREC1)*DSQRT(ER(IREC1,2))+DEERS(IREC1)), 2 XR(IREC1,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) 244 CONTINUE 242 CONTINUE C DO 344 IREC1=1,NREC2 DELR=DMAX1(1.0D-20,TS(IREC1)+DEER(IREC1)) TS(IREC1)=CVMGT(ER(IREC1,2)*TS(IREC1)/DELR,TS(IREC1) 1 ,DELR.GT.ER(IREC1,2)) DEER(IREC1)=CVMGT(ER(IREC1,2)*DEER(IREC1)/DELR,DEER(IREC1) 1 ,DELR.GT.ER(IREC1,2)) 344 CONTINUE C DO 252 IREC1=1,NREC2 C IF(XR(IREC1,2).LT.0.0.OR.XR(IREC1,2).GT.TT) GO TO 252 I=MAX0(MIN0(JFIX(X2(IREC1)/CW+1.D0),1000),1) C I=DMAX(MIN(INT(XR(IREC1,2)/CW+1.),100),1) DENTR(I)=DENTR(I)+TS(IREC1) DMGNR(I)=DMGNR(I)+T(IREC1) IONR(I)=IONR(I)+DEER(IREC1) C ELER(I,JJR(IREC1,2))=ELER(I,JJR(IREC1,2))+T(IREC1) C ELIR(I,JJR(IREC1,2))=ELIR(I,JJR(IREC1,2))+DEER(IREC1) IF(T(IREC1).LE.DI(JJR(IREC1,1))) GO TO 84 CC EPS(IV)=F1(JJJ(IV))*DEN(IV) CC G=EPS(IV)+.40244*EPS(IV)**.75+3.4008*EPS(IV)**.16667 CC MOT=DEN(IV)/(1.+K2(LLL(IV))*G) CC CASMOT(I)=CASMOT(I)+MOT ELGDR(I)=ELGDR(I)+T(IREC1) C ELDR(I,JJR(IREC1,2))=ELDR(I,JJR(IREC1,2))+T(IREC1) ICDR(I,JJR(IREC1,2))=ICDR(I,JJR(IREC1,2))+1 ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))= 1 ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))+1 GO TO 252 84 PHONR(I)=PHONR(I)+T(IREC1) C ELPR(I,JJJ(IREC1,2))=ELPR(I,JJJ(IREC1,2))+T(IREC1) 252 CONTINUE DO 82 IREC1=1,NREC2 ICDIR=ICDIR+JFIX(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1)))) ICSBR=ICSBR+JFIX(CVMGT(1.D0,0.D0,T(IREC1).GT.SB(1))) ICSUMR=ICSUMR+JFIX(CVMGT(1.D0,0.D0,TS(IREC1).GT.0.D0)) ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= 1 ICDIRJ(JJR(IREC1,2),JJR(IREC1,1)) 2 +JFIX(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1)))) 82 CONTINUE DO 253 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) 1 .OR.XR(IREC1,2).GT.SUT 2 .OR.(XR(IREC1,2).GT.RD.AND.XR(IREC1,2).LT.RT) 253 CONTINUE 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 cc IF(ER(IREC1,2).LE.SFE) GO TO 255 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(1001,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 C IF(SXR(I).EQ.0.) GO TO 182 C SXR(I)=DMAX1(SXR(I),1.0E-12) 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) C 182 CONTINUE 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 C MEAGSL(IESLOG,IGG,IAG)=MEAGSL(IESLOG,IGG,IAG)+1 C MEAGSL(IESLOG,IGG,21)=MEAGSL(IESLOG,IGG,21)+1 C MEAGSL(75,IGG,IAG)=MEAGSL(75,IGG,IAG)+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 C KADRIP(IAG,JJR(IREC1,2))=KADRIP(IAG,JJR(IREC1,2))+1 GO TO 88 86 ISPIST(JJR(IREC1,2))=ISPIST(JJR(IREC1,2))+1 ESPIST(JJR(IREC1,2))=ESPIST(JJR(IREC1,2))+ESPT C KADRIS(IAG,JJR(IREC1,2))=KADRIS(IAG,JJR(IREC1,2))+1 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 C KADROP(IAG,JJR(IREC1,2))=KADROP(IAG,JJR(IREC1,2))+1 GO TO 88 87 ISPOST(JJR(IREC1,2))=ISPOST(JJR(IREC1,2))+1 ESPOST(JJR(IREC1,2))=ESPOST(JJR(IREC1,2))+ESPT C KADROS(IAG,JJR(IREC1,2))=KADROS(IAG,JJR(IREC1,2))+1 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(1001,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 IF(ALPHA.LT.1.0) GO TO 181 C MEAGS(IESPT,IGG,IAGS) = MEAGS(IESPT,IGG,IAGS)+1 C MEAGS(102,IGG,IAGS) = MEAGS(102,IGG,IAGS)+1 C MEAGS(IESPT,IGG,22) = MEAGS(IESPT,IGG,22)+1 C MEAGSL(IESLOG,IGG,IAG)=MEAGSL(IESLOG,IGG,IAG)+1 C MEAGSL(IESLOG,IGG,21)=MEAGSL(IESLOG,IGG,21)+1 C MEAGSL(75,IGG,IAG)=MEAGSL(75,IGG,IAG)+1 C 181 CONTINUE 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 C 247 CONTINUE C 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 C IF(IH1.EQ.0.AND.IH.EQ.NH) GO TO 140 C C PROJECTILE CANDIDATE FOR REFLECTION C DO 29 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 29 CONTINUE 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 C IF(X(IV).LT.0..OR.X(IV).GT.TT) GO TO 110 IP = MAX0( MIN0( JFIX(PL(IV)/CW+1.D0), 1000), 1) IPL(IP)=IPL(IP)+1 I1 = MAX0( MIN0( JFIX(X(IV)/CW+1.D0), 1001), 0) IRP(I1)=IRP(I1)+1 c c Berechnung der gestoppten Teilchen im jeweiligen Layer c LowTiefe = 0.D0 UpTiefe = DX(1) c 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 c 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 C IJKLMN=IJKLMN+1 C ESVDL(IJKLMN)=ES 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( JFIX(PL(IV)/CW+1.D0), 1000), 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 C IF(ALPHA.LT.1.0) GO TO 183 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 C 183 CONTINUE 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( JFIX(PL(IV)/CW+1.D0), 1000), 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 C IF(ALPHA.LT.1.0) GO TO 183 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 C 183 CONTINUE 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 C 110 IF(IH.EQ.NH) GO TO 130 C 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)) CC703 CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,1) EMX = EMX+E(IV) CC WRITE(6,*) E(IV) GO TO 707 C 702 IF (EQUAL(Esig,0.D0)) THEN C FIXED PROJECTILE ENERGY C WRITE(*,*)' Da Esig=0 ist E=E0' 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 C WRITE(*,*)E(IV),Epar,E0 ENDIF C TAUPSI(IV)=0.D0 C C IF(ALPHA.EQ.-2.) GO TO 705 C IF(ALPHA.EQ.-1.) GO TO 706 IF(EQUAL(ALPHA,-2.D0)) GO TO 705 IF(EQUAL(ALPHA,-1.D0)) GO TO 706 C IF(EQUAL(ALPHASIG,0.D0))THEN C FIXED PROJECTILE ANGLE C WRITE(88,*)ALPHA,CALFA,SALFA 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) C WRITE(88,'(5(F12.5))')ALPHA,ALPHASIG,ALPHApar,CALFA,SALFA COSX(IV) = CALFA COSY(IV) = SALFA COSZ(IV) = 0.D0 SINE(IV) = COSY(IV) ENDIF C GO TO 707 C C COSINE ANGLE DISTRIBUTION C CC705 RPHI=PI2*RANF() CC705 RPHI=PI2*DRAND48() 705 RPHI=PI2*DBLE(RAN(ISEED)) CC RTHETA=RANF() CC RTHETA=DRAND48() RTHETA=DBLE(RAN(ISEED)) 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 C CC RPHI=PI2*RANF() CC RPHI=PI2*DRAND48() RPHI=PI2*DBLE(RAN(ISEED)) CC RTHETA=RANF() CC RTHETA=DRAND48() RTHETA=DBLE(RAN(ISEED)) 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 C CC709 RPHI=PI2*RANF() CC709 RPHI=PI2*DRAND48() 709 RPHI=PI2*DBLE(RAN(ISEED)) CC RTHETA=RANF() CC RTHETA=DRAND48() RTHETA=DBLE(RAN(ISEED)) 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 C 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) CC RA1=CVMGT(RANF(),1.,X0.LE.0.) CC RA1=CVMGT(DRAND48(),1.,X0.LE.0.) RA1=CVMGT(DBLE(RAN(ISEED)),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 CC IF(IVMIN.LE.1) GO TO 134 DO 128 IV=1,IVMIN-1 LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) 128 CONTINUE DO 129 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) 129 CONTINUE 134 CONTINUE C IF(IVMAX.LT.IVMIN) GO TO 132 CC IF(IVMAX.LT.IH1) GO TO 132 DO 133 IV=IVMAX+1,IH1 LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) 133 CONTINUE DO 131 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) 131 CONTINUE 132 CONTINUE C GO TO 1 C 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 c IF(ZT(1,2).LT.1.0D-3) THEN epsilon = 32.55D0*(MT(1,1)/M1)/(1.D0+(MT(1,1)/M1))* 1 1.D0/(Z1*ZT(1,1)*DSQRT(Z1**2.D0/3.D0+ZT(1,1)**2.D0/3.D0))* 2 E0keV prcoeff = prc(1)*DLOG(prc(2)*epsilon+DEXP(1.D0))/ 1 (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 CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3), 1 Date_Time) C IF(Date_Time(2).EQ.1) THEN month_stop='Jan.' days_total_stop=Date_Time(3) ELSEIF(Date_Time(2).EQ.2) THEN month_stop='Feb.' days_stop_total=31+Date_Time(3) ELSEIF(Date_Time(2).EQ.3) THEN month_stop='Mar.' days_stop_total=31+28+Date_Time(3) ELSEIF(Date_Time(2).EQ.4) THEN month_stop='Apr.' days_stop_total=31+28+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.5) THEN month_stop='May ' days_stop_total=31+28+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.6) THEN month_stop='Jun.' days_stop_total=31+28+31+30+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.7) THEN month_stop='Jul.' days_stop_total=31+28+31+30+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.8) THEN month_stop='Aug.' days_stop_total=31+28+31+30+31+30+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.9) THEN month_stop='Sep.' days_stop_total=31+28+31+30+31+30+31+31+Date_Time(3) ELSEIF(Date_Time(2).EQ.10) THEN month_stop='Oct.' days_stop_total=31+28+31+30+31+30+31+31+30+Date_Time(3) ELSEIF(Date_Time(2).EQ.11) THEN month_stop='Nov.' days_stop_total=31+28+31+30+31+30+31+31+30+31+Date_Time(3) ELSE month_stop='Dec.' days_stop_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3) ENDIF C READ(Real_Clock(1)(1:4),'(A4)')year_stop READ(Real_Clock(1)(7:8),'(A2)')day_stop READ(Real_Clock(2)(1:2),'(A2)')hour_stop READ(Real_Clock(2)(3:4),'(A2)')min_stop READ(Real_Clock(2)(5:6),'(A2)')sec_stop C C how many seconds are needed for the simulation ?? C seconds_stop_total=Date_Time(7)+(Date_Time(6)*60)+ 1 (Date_Time(5)*3600)+(days_stop_total-1)*86400 C WRITE(21,*) WRITE(21,10051)day_stop,month_stop,year_stop, 1 hour_stop,min_stop,sec_stop 10051 FORMAT(1x,' TrimSp simulation ended at: ',A2,'.',A4,1x,A4, 1 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') C 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 1 ,8HALPHASIG,7X,2HEF,7X 2 ,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 1 ,IPOTR,IRL 1406 FORMAT(/7X,2HNH,8X,2HRI,5X,3HRI2,5X,3HRI3,11X,2HX0,8X,2HRD,8X,2HCW 1 ,8X,2HCA 1 ,7X,3HKK0,3X,4HKK0R,3X,5HKDEE1,2X,5HKDEE2,2X,4HIPOT,3X,5HIPOTR 2 ,3X,3HIRL/I10,3F10.2,1F13.2,3F10.2,1X,7I7) WRITE(21,1408) 1408 FORMAT(//13X,2HDX,6X,3HRHO,4X,2HCK,2X 1 ,5HZ(,1),1X,5HZ(,2),1X,5HZ(,3),1X,5HZ(,4),1X,5HZ(,5),2X 2 ,5HM(,1),2X,5HM(,2),2X,5HM(,3),2X,5HM(,4),2X,5HM(,5),1X 3 ,5HC(,1),1X,5HC(,2),1X,5HC(,3),1X,5HC(,4),1X,5HC(,5)) DO 1410 I=1,L WRITE(21,1412) I,DX(I),RHO(I),CK(I),(ZT(I,J),J=1,5) 1 ,(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) 1410 CONTINUE WRITE(21,1414) 1414 FORMAT(//27X,'***',2X,'SBE(LAYER,ELEMENT)',2X,'***',5X 1 ,'***',5X,'ED(LAYER,ELEMENT)',5X,'***',5X 2 ,'***',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,' 1 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,' 1 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 1 ,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 1 ,'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) WRITE(22,14280)rgenam 14280 FORMAT(1H1,//30X,'* RANGE DATA *',5X,A12) 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 1 ,'SFE',6X,'INEL',9X,'L',8X,'LJ'/ 2 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 1 ,'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)' 1 ,7X,'LM(I)',5X,'PDMAX(I)',5X,'ASIG(I)',7X,'SB(I)',7X,'XX(I)' 2 ,8X,'NJ(I)') DO 472 I=1,L WRITE(21,473) I,EPS0(I),Z2(I),M2(I),ARHO(I),LM(I),PDMAX(I),ASIG(I) 1 ,SB(I),XX(I),NJ(I) 473 FORMAT(/1X,I1,6H.LAYER,1X,9E12.4,I10) 472 CONTINUE WRITE(21,474) 474 FORMAT(//13X, 1 'A1(1)',3X,'A1(2)',3X,'A1(3)',3X,'A1(4)',3X,'A1(5)',3X, 1 'A1(6)',3X,'A1(7)',3X,'A1(8)',3X,'A1(9)',2X,'A1(10)',2X, 1 'A1(11)',2X,'A1(12)',2X,'A1(13)',2X,'A1(14)',2X,'A1(15)',2X, 1 'A1(16)',2X,'A1(17)',2X,'A1(18)',2X,'A1(19)',2X,'A1(20)',2X, 1 'A1(21)',2X,'A1(22)',2X,'A1(23)',2X,'A1(24)',2X,'A1(25)',2X, 1 'A1(26)',2X,'A1(27)',2X,'A1(28)',2X,'A1(29)',2X,'A1(30)',2X, 1 '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, 1 'KOR1(1)',1X,'KOR1(2)',1X,'KOR1(3)',1X,'KOR1(4)',1X,'KOR1(5)', 2 1X,'KOR1(6)',1X,'KOR1(7)',1X,'KOR1(8)',1X,'KOR1(9)',1X,'KOR1(A)', 3 1X,'KOR1(B)',1X,'KOR1(C)',1X,'KOR1(D)',1X,'KOR1(E)',1X,'KOR1(F)', 4 1X,'KOR1(G)',1X,'KOR1(H)',1X,'KOR1(I)',1X,'KOR1(J)',1X,'KOR1(K)', 5 1X,'KOR1(L)',1X,'KOR1(M)',1X,'KOR1(N)',1X,'KOR1(O)',1X,'KOR1(P)', 6 1X,'KOR1(Q)',1X,'KOR1(R)',1X,'KOR1(S)',1X,'KOR1(T)',1X,'KOR1(U)', 7 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, 1 'A(I,1)',2X,'A(I,2)',2X,'A(I,3)',2X,'A(I,4)',2X,'A(I,5)',2X, 2 'A(I,6)',2X,'A(I,7)',2X,'A(I,8)',2X,'A(I,9)',1X,'A(I,10)',1X, 3 'A(I,11)',1X,'A(I,12)',1X,'A(I,13)',1X,'A(I,14)',1X,'A(I,15)',1X, 4 'A(I,16)',1X,'A(I,17)',1X,'A(I,18)',1X,'A(I,19)',1X,'A(I,20)',1X, 5 'A(I,21)',1X,'A(I,22)',1X,'A(I,23)',1X,'A(I,24)',1X,'A(I,25)',1X, 6 'A(I,26)',1X,'A(I,27)',1X,'A(I,28)',1X,'A(I,29)',1X,'A(I,30)',1X, 7 'A(I,31)',1X,'A(I,32)',1X,'A(I,33)',1X,'A(I,34)',1X,'A(I,35)') DO 478 I=1,LJ WRITE(21,477) (A(I,J),J=1,LJ) 477 FORMAT(/1X,9X,35F8.5) 478 CONTINUE WRITE(21,490) 490 FORMAT(//11X, 1 'KOR(,1)',1X,'KOR(,2)',1X,'KOR(,3)',1X,'KOR(,4)',1X,'KOR(,5)',1X, 2 'KOR(,6)',1X,'KOR(,7)',1X,'KOR(,8)',1X,'KOR(,9)',1X,'KOR(,A)',1X, 3 'KOR(,B)',1X,'KOR(,C)',1X,'KOR(,D)',1X,'KOR(,E)',1X,'KOR(,F)',1X, 4 'KOR(,G)',1X,'KOR(,H)',1X,'KOR(,I)',1X,'KOR(,J)',1X,'KOR(,K)',1X, 5 'KOR(,L)',1X,'KOR(,M)',1X,'KOR(,N)',1X,'KOR(,O)',1X,'KOR(,P)',1X, 6 'KOR(,Q)',1X,'KOR(,R)',1X,'KOR(,S)',1X,'KOR(,T)',1X,'KOR(,U)',1X, 7 'KOR(,V)',1X,'KOR(,W)',1X,'KOR(,X)',1X,'KOR(,Y)',1X,'KOR(,Z)') DO 491 I=1,LJ WRITE(21,492) (KOR(I,J),J=1,LJ) 492 FORMAT(/1X,9X,35F8.5) 491 CONTINUE C WRITE(6,479) C 479 FORMAT(//13X, C 1 'F(I,1)',6X,'F(I,2)',6X,'F(I,3)',6X,'F(I,4)',6X,'F(I,5)',5X, C 2 'KOR(I,1)',4X,'KOR(I,2)',4X,'KOR(I,3)',4X,'KOR(I,4)',4X, C 3 'KOR(I,5)') C DO 480 I=1,L C WRITE(6,481) I,(F(I,J),J=1,5),(KOR(I,J),J=1,5) C 481 FORMAT(/1X,I1,6H.LAYER,1X,10E12.4) C 480 CONTINUE C WRITE(6,483) C 483 FORMAT(//12X, C 1 'KL(I,1)',5X,'KL(I,2)',5X,'KL(I,3)',5X,'KL(I,4)',5X,'KL(I,5)',6X, C 2 'K(I,1)',6X,'K(I,2)',6X,'K(I,3)',6X,'K(I,4)',6X,'K(I,5)') C DO 482 I=1,L C WRITE(6,485) I,(KL(I,J),J=1,5),(K(I,J),J=1,5) C 485 FORMAT(/1X,I1,6H.LAYER,1X,10E12.4) C 482 CONTINUE C C INTEGRAL IMPLANTATION , SPUTTERING , BACKSCATTERING , TRANSMISSION C IIM=NH-IB-IT YH=DFLOTJ(IIM) HN=DFLOTJ(NH) EMV=CVMGT(EMX/HN,E0,E0.LE.0.D0) EIM=DFLOTJ(IIM)*EMV DO 1550 J=1,LJ ISPA = ISPA+IBSP(J) 1550 ESPA = ESPA+EBSP(J) DO 1702 J=1,LJ ISPAT = ISPAT+ITSP(J) 1702 ESPAT = ESPAT+ETSP(J) WRITE(21,1500) IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,ISPAT,ESPAT, 1 tryE,negE,epsilon,prcoeff 1500 FORMAT(1H1,//11X,20HIMPLANTED PARTICLES=,I7,5X,7HENERGY=,E10.4, 1 3H EV/7X,24HBACKSCATTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, 2 3H EV/9X,22HTRANSMITTED PARTICLES=,I7,5X,7HENERGY=,E10.4, 3 3H EV/7X,24HBACKSPUTTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, 4 3H EV/6X,'TRANSM. SPUTT. PARTICLES=',I7,5X,7HENERGY=,E10.4, 5 3H EV/15X,16HTRIED PARTICLES=,I7 6 /9X,22HPARTICLES with neg. E=,I7, 7 /6X,25HTHOMAS FERMI RED. ENERGY=,2X,E10.4, 8 /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=DFLOTJ(ICDI)/HN CST=DFLOTJ(ICSUM-ICDI) WRITE(21,1511) AVCSUM,AVCDIS,AVCSMS 1511 FORMAT(//2X,'PROJECTILES : ', 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F8.3/65X, 3 '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, 1 XSUM,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,YH) CALL MOMENTS(FIR0,SER,THR,FOR,FIR,SIR,SIGMAR,DFIR0,DSER,DTHR, 1 RSUM,R2SUM,R3SUM,R4SUM,R5SUM,R6SUM,YH) CALL MOMENTS(FIP0,SEP,THP,FOP,FIP,SIP,SIGMAP,DFIP0,DSEP,DTHP, 1 PLSUM,PL2SUM,PL3SUM,PL4SUM,PL5SUM,PL6SUM,YH) CALL MOMENTS(FIE0,SEE,THE,FOE,FIE,SIE,SIGMAE,DFIE0,DSEE,DTHE, 1 EEL,EEL2,EEL3,EEL4,EEL5,EEL6,CSUM) CALL MOMENTS(FIW0,SEW,THW,FOW,FIW,SIW,SIGMAW,DFIW0,DSEW,DTHW, 1 EELWC,EELWC2,EELWC3,EELWC4,EELWC5,EELWC6,CSUM) CALL MOMENTS(FII0,SEI,THI,FOI,FII,SII,SIGMAI,DFII0,DSEI,DTHI, 1 EIL,EIL2,EIL3,EIL4,EIL5,EIL6,CSUM) CALL MOMENTS(FIS0,SES,THS,FOS,FIS,SIS,SIGMAS,DFIS0,DSES,DTHS, 1 EPL,EPL2,EPL3,EPL4,EPL5,EPL6,CST) WRITE(21,7117) 7117 FORMAT(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',4X, 1 ' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',3X,' ERROR 2.M ', 2 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 c IF(YH.LT.1.D0) GO TO 7235 CALL MOMENT(X1SD,X2SD,X3SD,X4SD,X5SD,X6SD 1 ,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=DFLOTJ(ICSUMR) ACSUMR=CSUMR/HN ACDISR=DFLOTJ(ICDIR)/HN ACSBER=DFLOTJ(ICSBR)/HN WRITE(21,1599) ACSUMR,ACDISR,ACSBER 1599 FORMAT(///2X,'RECOILES (PER PROJ.) : ', 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/76X, 3 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3) IF(NPA+NSA.EQ.0) GO TO 1453 ACSUR=CSUMR/(DFLOTJ(NPA+NSA)) ACDIR=DFLOTJ(ICDIR)/(NPA+NSA) ACSBR=DFLOTJ(ICSBR)/(NPA+NSA) WRITE(21,1598) ACSUR,ACDIR,ACSBR 1598 FORMAT(/2X,'RECOILES (PER KNOCKON) : ', 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.3,3X, 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/,76X, 3 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3/) IF(NJ(1).LT.2) GO TO 1453 ACDR11=DFLOTJ(ICDIRJ(1,1))/(NPA+NSA) ACDR12=DFLOTJ(ICDIRJ(1,2))/(NPA+NSA) ACDR21=DFLOTJ(ICDIRJ(2,1))/(NPA+NSA) ACDR22=DFLOTJ(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/ 1 ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-2 = ',F10.3/ 2 ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-1 = ',F10.3/ 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, 110HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X,10HELLOSS(EV), 22X,10HDAMAGE(EV),2X,10HPHONON(EV),2X,10HCASCAD(EV),5X,3HDPA/) C IF(depth_interval_flag.EQ.0) THEN WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT AGREE WIT 1H LAYER THICKNESS' WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT AGREE WIT 1H LAYER THICKNESS' ENDIF C WRITE(22,6002) 6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES) IF(YH.LT.1.D0) GO TO 603 DO 602 I=0,1001 RIRP(I) = DFLOTJ(IRP(I))/YH 602 CONTINUE 603 D1=0. D2=CW WRITE(21,601) D1,IRP(0),RIRP(0) 601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4) c DO 1441 J=1,NJ(1) DO 1441 J=1,LJ DO 1441 I=1,1000 ICDT(I)=ICDT(I)+ICD(I,J) ICDTR(I)=ICDTR(I)+ICDR(I,J) 1441 CONTINUE DO 1442 K=1,NJ(1) DO 1442 J=1,NJ(1) DO 1442 I=1,1000 ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) 1442 CONTINUE DO 35 I=0,1001 IIRP=IIRP+IRP(I) TRIRP=TRIRP+RIRP(I) 35 CONTINUE DO 36 I=1,1000 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) C TCASMOR=TCASMOR+CASMOTR(I) ICDTTR=ICDTTR+ICDTR(I) 36 CONTINUE do im1=1000,1,-1 C if(ipl(im1).ne.0.or.ion(im1).ne.0.) go to 20 if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 enddo im1=1 20 im1=min0(im1+2,1000) DO 11 I=1,im1 WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I),DMGN(I) 1 ,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(1001),RIRP(1001) 604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4) WRITE(21,710) IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO 1 ,ICDTT 710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8) DO 1531 J=1,NJ(1) DO 1531 I=1,1000 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) C ELETR(J)=ELETR(J)+ELE(I,J) C ELITR(J)=ELITR(J)+ELI(I,J) C ELPTR(J)=ELPTR(J)+ELP(I,J) C ELDTR(J)=ELDTR(J)+ELD(I,J) 1531 CONTINUE c IF(NJ(1).LT.2) GO TO 1455 WRITE(21,1521) 1521 FORMAT(1H1,4X,'DEPTH(A)' 1 ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)' 2 ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)' 3 ,2X,'DPA(1)',2X,'DPA(2)'/) D1=0. D2=CW do im2=1000,1,-1 C if(eli(im2,1).ne.0..or.eli(im2,2).ne.0.) go to 30 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,1000) DO 1525 I=1,im2 WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1) 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) 1 ,ELIT(2),ELET(2),ELDT(2),ELPT(2),ICDJT(1),ICDJT(2) 1533 FORMAT(/14X,1P8E12.4,2I8///) C WRITE(21,1481) C1481 FORMAT(1H1,2X,'D(A)' C 1 ,2X,' DAMAGE(1)',3X,' DAMAGE(2)',1X,' DAMAGE(3)' C 2 ,2X,' DAMAGE(4)',3X,' DAMAGE(5)'/) C DO 1482 I=1,100 C WRITE(6,1483) I,ELD(I,1),ELD(I,2),ELD(I,3),ELD(I,4),ELD(I,5) C1482 CONTINUE C1483 FORMAT(1X,I5,1P5E14.4) C WRITE(6,1484) ELDT(1),ELDT(2),ELDT(3),ELDT(4),ELDT(5) C1484 FORMAT(/1X,5X,1P5E14.4///) DO 1491 I=1,L-1 ILD(I)=IDINT(XX(I)/CW+0.01D0) IF(ILD(I).GT.1000) ILD(I)=1000 DO 1492 J=1,ILD(I) DLI(I)=DLI(I)+DMGN(J) 1492 CONTINUE 1491 CONTINUE DLI(L)=TDMGN C WRITE(21,*) 'L=',L,' XX=',XX,' DLI=',DLI 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), 1 5X,10HINLOSS(EV),3X,10HTELOSS(EV),3X,10HELLOSS(EV), 2 3X,10HDAMAGE(EV),3X,10HPHONON(EV),5X,3HDPA, c 3 2X,6HDPA(1),2X,6HDPA(2)/) 3 2X,6HDPA(1),2X,6HDPA(2), 4 1X,5H(1-1),1X,5H(1-2),1X,5H(2-1),1X,5H(2-2)/) D1=0.D0 D2=CW do im3=1000,1,-1 if (.not.equal(ionr(im3),0.D0)) go to 31 C if(ionr(im3).ne.0.) goto 31 enddo im3=1 31 im3=MIN0(im3+2,1000) DO 1594 I=1,im3 WRITE(21,1595) D1,D2,IONR(I),DENTR(I),DMGNR(I),ELGDR(I),PHONR(I) 1 ,ICDTR(I),ICDIRN(I,1),ICDIRN(I,2) 2 ,ICDIRI(I,1,1),ICDIRI(I,1,2),ICDIRI(I,2,1),ICDIRI(I,2,2) c1595 FORMAT(1X,F6.0,1H-,F6.0,1P1E14.4,4E13.4,3I8) 1595 FORMAT(1X,F6.0,1H-,F6.0,1P1E14.4,4E13.4,3I8,4I6) D1=D2 1594 D2=D2+CW WRITE(21,1596) TIONR,TDENTR,TDMGNR,TELGDR,TPHONR 1,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=DFLOTJ(IB) BIL=DFLOTJ(IBL) RN=BI/HN RE=EB/(HN*EMV) EMEANR=RE/RN EMEAN=EB/BI AVEB=EMEAN IF (equal(BI,1.0d0))GO TO 1506 C IF(BI.EQ.1.) 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.=' 1 ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' 2 ,1E11.4) IF(IB.EQ.0) GO TO 1512 CALL MOMENT(EB1B,EB2B,EB3B,EB4B,EB5B,EB6B 1 ,EB,EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM,BI) CALL MOMENT(EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL 1 ,EB1SUL,EB2SUL,EB3SUL,EB4SUL,EB5SUL,EB6SUL,BIL) CALL MOMENTS(FIB0,SEB,THB,FOB,FIB,SIB,SIGMAB,DFIB0,DSEB,DTHB, 1 EB,EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM,BI) CALL MOMENT(PL1S,PL2S,PL3S,PL4S,PL5S,PL6S 1 ,PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB,BI) CALL MOMENTS(FIPB0,SEPB,THPB,FOPB,FIPB,SIPB,SIGMPB 1 ,DFIPB0,DSEPB,DTHPB, 2 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*DFLOTJ(I) IF(IB.EQ.0) GO TO 1512 WRITE(21,1514) 1514 FORMAT(//5X,'POLAR ANGULAR DISTRIBUTION OF BACKSCATTERED PROJECTIL 1ES'//) DO 1516 I=1,20 1516 RKADB(I)=DFLOTJ(KADB(I))*20.D0/DFLOTJ(IB) 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=DFLOTJ(IT) TN=TIT/HN TE=ET/(HN*E0) TMEANR=TE/TN EMEANT=TMEANR*E0 IF (equal(TIT,1.0D0)) GO TO 1520 C IF(TIT.EQ.1.) 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.=' 1 ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' 2 ,1E11.4) CALL MOMENTS(FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT, 1 ET,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM,TIT) CALL MOMENTS(FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT 1 ,DFIPT0,DSEPT,DTHPT, 2 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' 1//) DO 1528 I=1,20 1528 RKADT(I)=DFLOTJ(KADT(I))*20.D0/DFLOTJ(IT) 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 1552 J=1,NJ(1) ISPAL(1) = ISPAL(1)+IBSP(J) 1552 ESPAL(1) = ESPAL(1)+EBSP(J) DO 1554 J=NJ(1)+1,JT(3) ISPAL(2) = ISPAL(2)+IBSP(J) 1554 ESPAL(2) = ESPAL(2)+EBSP(J) DO 1556 J=JT(3)+1,LJ ISPAL(3) = ISPAL(3)+IBSP(J) 1556 ESPAL(3) = ESPAL(3)+EBSP(J) WRITE(21,1558) ISPA,ESPA 1558 FORMAT(///,8X,'ALL SPUTTERED PARTICLES = ',I7,3X 1 ,'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 1 ,'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 1 ,'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV') 1562 CONTINUE IF(ISPA.EQ.0) GO TO 1700 DO 1572 J=1,LJ RIP(J)=DFLOTJ(ISPIP(J))/DFLOTJ(ISPA) RIS(J)=DFLOTJ(ISPIS(J))/DFLOTJ(ISPA) ROP(J)=DFLOTJ(ISPOP(J))/DFLOTJ(ISPA) ROS(J)=DFLOTJ(ISPOS(J))/DFLOTJ(ISPA) REIP(J)=ESPIP(J)/ESPA REIS(J)=ESPIS(J)/ESPA REOP(J)=ESPOP(J)/ESPA REOS(J)=ESPOS(J)/ESPA 1572 CONTINUE DO 1584 J=1,LJ IF(IBSP(J).EQ.0) GO TO 1584 RIPJ(J)=DFLOTJ(ISPIP(J))/DFLOTJ(IBSP(J)) RISJ(J)=DFLOTJ(ISPIS(J))/DFLOTJ(IBSP(J)) ROPJ(J)=DFLOTJ(ISPOP(J))/DFLOTJ(IBSP(J)) ROSJ(J)=DFLOTJ(ISPOS(J))/DFLOTJ(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)/DFLOTJ(ISPIP(J)) 3571 IF(ISPIS(J).EQ.0) GO TO 3572 ESPMIS(J)=ESPIS(J)/DFLOTJ(ISPIS(J)) 3572 IF(ISPOP(J).EQ.0) GO TO 3573 ESPMOP(J)=ESPOP(J)/DFLOTJ(ISPOP(J)) 3573 IF(ISPOS(J).EQ.0) GO TO 1571 ESPMOS(J)=ESPOS(J)/DFLOTJ(ISPOS(J)) 1571 CONTINUE 1573 CONTINUE DO 1578 J=1,LJ SPY(J)=DFLOTJ(IBSP(J))/HN 1578 SPE(J)=EBSP(J)/(HN*EMV) DO 1579 J=1,LJ IF (equal(SPY(J),0.0D0))GO TO 1579 C IF(SPY(J).EQ.0.0) 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 1574 J=1,NJ(1) WRITE(21,1576) J,ISPIP(J),RIP(J),RIPJ(J),ESPIP(J),REIP(J),REIPJ(J) 1 ,ESPMIP(J) 2 ,J,ISPIS(J),RIS(J),RISJ(J),ESPIS(J),REIS(J),REISJ(J) 3 ,ESPMIS(J) 4 ,J,ISPOP(J),ROP(J),ROPJ(J),ESPOP(J),REOP(J),REOPJ(J) 5 ,ESPMOP(J) 6 ,J,ISPOS(J),ROS(J),ROSJ(J),ESPOS(J),REOS(J),REOSJ(J) 7 ,ESPMOS(J) 1576 FORMAT(/9X,'ION IN , PRIMARY KO(',I1,') = ',I7,2F9.4,4X 1 ,'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/ 2 9X,'ION IN , SECOND. KO(',I1,') = ',I7,2F9.4,4X 3 ,'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/ 4 8X,'ION OUT , PRIMARY KO(',I1,') = ',I7,2F9.4,4X 5 ,'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4/ 6 8X,'ION OUT , SECOND. KO(',I1,') = ',I7,2F9.4,4X 7 ,'ENERGY = ',E10.4,' EV',2F9.4,4X,'MEAN ENERGY = ',E10.4) 1574 CONTINUE 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, 1 ' SPUTTERED ENERGY(',I1,') = ',1E10.3, 2 ' REL.MEAN ENERGY(',I1,') = ',1E10.3, 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 1 ,DFIES0,DSEES,DTHES, 2 EBSP1,EBSP2,EBSP3,EBSP4,EBSP5,EBSP6 3 ,EBSP(J),SPE2S(J),SPE3S(J),SPE4S(J),SPE5S(J) 4 ,SPE6S(J),YSP) CALL MOMENTN(FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,SIGMSL 1 ,DFIESL,DSEESL,DTHESL, 2 EBSP1L,EBSP2L,EBSP3L,EBSP4L,EBSP5L,EBSP6L 3 ,SPE1SL(J),SPE2SL(J),SPE3SL(J),SPE4SL(J),SPE5SL(J) 4 ,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 1 ,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 ' 1 ,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 1 ,'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) 1 ,REIPJ(J),ESPMIP(J) 2 ,J-NJ(1),ISPIS(J),RIS(J),RISJ(J),ESPIS(J),REIS(J) 3 ,REISJ(J),ESPMIS(J) 4 ,J-NJ(1),ISPOP(J),ROP(J),ROPJ(J),ESPOP(J),REOP(J) 5 ,REOPJ(J),ESPMOP(J) 6 ,J-NJ(1),ISPOS(J),ROS(J),ROSJ(J),ESPOS(J),REOS(J) 7 ,REOSJ(J),ESPMOS(J) 1586 CONTINUE WRITE(21,1577) DO 1592 J=NJ(1)+1,JT(3) WRITE(21,1582) J-NJ(1),SPY(J),J,SPE(J),J,REY(J),J,EMSP(J) 1592 CONTINUE 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 1 ,DFIES0,DSEES,DTHES, 2 EBSP(J),SPE2S(J),SPE3S(J),SPE4S(J),SPE5S(J) 3 ,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 1 ,'SPUTTERED ENERGY(',I1,') = ',E10.4,' EV') 1536 CONTINUE DO 1540 J=JT(3)+1,LJ WRITE(21,1576) J-JT(3),ISPIP(J),RIP(J),RIPJ(J),ESPIP(J),REIP(J) 1 ,REIPJ(J),ESPMIP(J) 2 ,J-JT(3),ISPIS(J),RIS(J),RISJ(J),ESPIS(J),REIS(J) 3 ,REISJ(J),ESPMIS(J) 4 ,J-JT(3),ISPOP(J),ROP(J),ROPJ(J),ESPOP(J),REOP(J) 5 ,REOPJ(J),ESPMOP(J) 6 ,J-JT(3),ISPOS(J),ROS(J),ROSJ(J),ESPOS(J),REOS(J) 7 ,REOSJ(J),ESPMOS(J) 1540 CONTINUE WRITE(21,1577) DO 1542 J=JT(3)+1,LJ WRITE(21,1582) J-JT(3),SPY(J),J-JT(3),SPE(J),J-JT(3),REY(J) 1 ,J-JT(3),EMSP(J) 1542 CONTINUE 1532 CONTINUE C C BACKWARD SPUTTERING : ANGULAR DISTRIBUTIONS C WRITE(21,1601) 1601 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF ALL BACKWARD SPUTTERED 1 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 1602 I=1,20 DO 1602 J=1,NJ(1) 1602 KADSL(I,1)=KADSL(I,1)+KADSJ(I,J) DO 1604 I=1,20 DO 1604 J=NJ(1)+1,JT(3) 1604 KADSL(I,2)=KADSL(I,2)+KADSJ(I,J) 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 ; 1LAYER 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) 1 ,(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 ; 1LAYER 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) 1 ,(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 ; 1LAYER 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) 1 ,(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 ; 1LAYER 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) 1 ,(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 1706 J=1,NJ(1) ISPALT(1) = ISPALT(1)+ITSP(J) 1706 ESPALT(1) = ESPALT(1)+ETSP(J) DO 1708 J=NJ(1)+1,JT(3) ISPALT(2) = ISPALT(2)+ITSP(J) 1708 ESPALT(2) = ESPALT(2)+ETSP(J) DO 1710 J=JT(3)+1,LJ ISPALT(3) = ISPALT(3)+ITSP(J) 1710 ESPALT(3) = ESPALT(3)+ETSP(J) WRITE(21,1712) ISPAT,ESPAT 1712 FORMAT(///,8X,'ALL SPUTTERED PARTICLES = ',I7,3X 1 ,'TOTAL SPUTTERED ENERGY = ',E10.4,3H EV//) DO 1711 J=1,L WRITE(21,1713) J,ISPALT(J),ESPALT(J) 1713 FORMAT(8X,'SPUTTERED PARTICLES (LAYER ',I1,') = ',I7,3X 1 ,'SPUTTERED ENERGY = ',E10.4,3H EV) 1711 CONTINUE DO 1732 J=1,LJ RIPT(J)=DFLOTJ(ISPIPT(J))/DFLOTJ(ISPAT) RIST(J)=DFLOTJ(ISPIST(J))/DFLOTJ(ISPAT) ROPT(J)=DFLOTJ(ISPOPT(J))/DFLOTJ(ISPAT) ROST(J)=DFLOTJ(ISPOST(J))/DFLOTJ(ISPAT) REIPT(J)=ESPIPT(J)/ESPAT REIST(J)=ESPIST(J)/ESPAT REOPT(J)=ESPOPT(J)/ESPAT 1732 REOST(J)=ESPOST(J)/ESPAT 1715 CONTINUE DO 1717 J=1,LJ IF(ISPIPT(J).EQ.0) GO TO 4571 ESPMIPT(J)=ESPIPT(J)/DFLOTJ(ISPIPT(J)) 4571 IF(ISPIST(J).EQ.0) GO TO 4572 ESPMIST(J)=ESPIST(J)/DFLOTJ(ISPIST(J)) 4572 IF(ISPOPT(J).EQ.0) GO TO 4573 ESPMOPT(J)=ESPOPT(J)/DFLOTJ(ISPOPT(J)) 4573 IF(ISPOST(J).EQ.0) GO TO 1717 ESPMOST(J)=ESPOST(J)/DFLOTJ(ISPOST(J)) 1717 CONTINUE DO 1736 J=1,LJ SPYT(J)=DFLOTJ(ITSP(J))/DFLOTJ(NH) 1736 SPET(J)=ETSP(J)/(NH*E0) DO 1737 J=1,LJ IF (equal(SPYT(J),0.0D0))GO TO 1737 C IF(SPYT(J).EQ.0.0) 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 1716 J=1,NJ(1) WRITE(21,1564) J,ITSP(J),J,ETSP(J) 1716 CONTINUE DO 1734 J=1,NJ(1) WRITE(21,1581) J,ISPIPT(J),RIPT(J),ESPIPT(J),REIPT(J),ESPMIPT(J) 1 ,J,ISPIST(J),RIST(J),ESPIST(J),REIST(J),ESPMIST(J) 2 ,J,ISPOPT(J),ROPT(J),ESPOPT(J),REOPT(J),ESPMOPT(J) 3 ,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 1 ,'ENERGY = ',E10.4,' EV',1F9.4,4X,'MEAN ENERGY = ',E10.4/ 2 9X,'ION IN , SECOND. KO(',I1,') = ',I7,1F9.4,4X 3 ,'ENERGY = ',E10.4,' EV',1F9.4,4X,'MEAN ENERGY = ',E10.4/ 4 8X,'ION OUT , PRIMARY KO(',I1,') = ',I7,1F9.4,4X 5 ,'ENERGY = ',E10.4,' EV',1F9.4,4X,'MEAN ENERGY = ',E10.4/ 6 8X,'ION OUT , SECOND. KO(',I1,') = ',I7,1F9.4,4X 7 ,'ENERGY = ',E10.4,' EV',1F9.4,4X,'MEAN ENERGY = ',E10.4) WRITE(21,1577) DO 1738 J=1,NJ(1) WRITE(21,1582) J,SPYT(J),J,SPET(J),J,REYT(J),J,EMSPT(J) 1738 CONTINUE 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 1722 J=NJ(1)+1,JT(3) WRITE(21,1564) J-NJ(1),ITSP(J),J-NJ(1),ETSP(J) 1722 CONTINUE DO 1746 J=NJ(1)+1,JT(3) WRITE(21,1581) J-NJ(1),ISPIPT(J),RIPT(J),ESPIPT(J),REIPT(J) 1 ,ESPMIPT(J) 2 ,J-NJ(1),ISPIST(J),RIST(J),ESPIST(J),REIST(J) 3 ,ESPMIST(J) 4 ,J-NJ(1),ISPOPT(J),ROPT(J),ESPOPT(J),REOPT(J) 5 ,ESPMOPT(J) 6 ,J-NJ(1),ISPOST(J),ROST(J),ESPOST(J),REOST(J) 7 ,ESPMOST(J) 1746 CONTINUE WRITE(21,1577) DO 1748 J=NJ(1)+1,JT(3) WRITE(21,1582) J-NJ(1),SPYT(J),J-NJ(1),SPET(J),J-NJ(1),REYT(J) 1 ,J-NJ(1),EMSPT(J) 1748 CONTINUE 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 1728 J=JT(3)+1,LJ WRITE(21,1564) J-JT(3),ITSP(J),J-JT(3),ETSP(J) 1728 CONTINUE DO 1750 J=JT(3)+1,LJ WRITE(21,1581) J-JT(3),ISPIPT(J),RIPT(J),ESPIPT(J),REIPT(J) 1 ,ESPMIPT(J) 2 ,J-JT(3),ISPIST(J),RIST(J),ESPIST(J),REIST(J) 3 ,ESPMIST(J) 4 ,J-JT(3),ISPOPT(J),ROPT(J),ESPOPT(J),REOPT(J) 5 ,ESPMOPT(J) 6 ,J-JT(3),ISPOST(J),ROST(J),ESPOST(J),REOST(J) 7 ,ESPMOST(J) 1750 CONTINUE WRITE(21,1577) DO 1752 J=JT(3)+1,LJ WRITE(21,1582) J-JT(3),SPYT(J),J-JT(3),SPET(J),J-JT(3),REYT(J) 1 ,J-JT(3),EMSPT(J) 1752 CONTINUE 1749 CONTINUE C C TRANSMISSION SPUTTERING : ANGULAR DISTRIBUTIONS C WRITE(21,1760) 1760 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF ALL TRANSMISSION SPUTT 1ERED PARTICLES'//) DO 1762 I=1,20 1762 RKADST(I)=KADST(I)*20.D0/ISPAT 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 1766 I=1,20 DO 1768 J=1,NJ(1) 1768 KDSTL(I,1)=KDSTL(I,1)+KDSTJ(I,J) DO 1770 J=NJ(1)+1,JT(3) 1770 KDSTL(I,2)=KDSTL(I,2)+KDSTJ(I,J) 1766 CONTINUE DO 1753 J=1,2 IF(ISPAL(J).EQ.0) GO TO 1754 DO 1772 I=1,20 1772 RKDSTL(I,J)=KDSTL(I,J)*20.D0/ISPAL(J) 1754 CONTINUE 1753 CONTINUE DO 1755 J=1,JT(3) IF(ITSP(J).EQ.0) GO TO 1756 DO 1774 I=1,20 1774 RKDSTJ(I,J)=KDSTJ(I,J)*20.D0/ITSP(J) 1756 CONTINUE 1755 CONTINUE WRITE(21,1776) 1776 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED PARTICLES ; 1LAYER 1'//) WRITE(21,1518) (AI(I),I=1,20),(KDSTL(I,1),I=1,20) 1 ,(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 ; 1LAYER 1 , SPECIES ',I1//) WRITE(21,1518) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) 1 ,(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 ; 1LAYER 2'//) WRITE(21,1518) (AI(I),I=1,20),(KDSTL(I,2),I=1,20) 1 ,(RKDSTL(I,2),I=1,20) IF(NJ(2).EQ.1) GO TO 1800 DO 1788 J=NJ(1)+1,JT(3) WRITE(21,1790) J-NJ(1) 1790 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED PARTICLES ; 1LAYER 2 , SPECIES ',I1//) WRITE(21,1518) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) 1 ,(RKDSTJ(I,J),I=1,20) 1788 CONTINUE GO TO 1800 1764 DO 1761 I=1,20 DO 1763 J=1,NJ(2) 1763 KDSTL(I,1)=KDSTL(I,1)+KDSTJ(I,J) DO 1765 J=NJ(2)+1,LJ-NJ(1) 1765 KDSTL(I,2)=KDSTL(I,2)+KDSTJ(I,J) 1761 CONTINUE DO 1799 J=1,2 IF(ISPALT(J+1).EQ.0) GO TO 1799 DO 1767 I=1,20 1767 RKDSTL(I,J)=KDSTL(I,J)*20.D0/ISPALT(J+1) 1799 CONTINUE DO 1797 J=1,LJ-NJ(1) IF(ITSP(J+NJ(1)).EQ.0) GO TO 1797 DO 1769 I=1,20 1769 RKDSTJ(I,J)=KDSTJ(I,J)*20.D0/ITSP(J+NJ(1)) 1797 CONTINUE WRITE(21,1771) 1771 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED PARTICLES ; 1LAYER 2'//) WRITE(21,1518) (AI(I),I=1,20),(KDSTL(I,1),I=1,20) 1 ,(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 ; 1LAYER 2 ; SPECIES ',I1//) WRITE(21,1518) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) 1 ,(RKDSTJ(I,J),I=1,20) 1775 CONTINUE 1773 CONTINUE WRITE(21,1779) 1779 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED PARTICLES ; 1LAYER 3'//) WRITE(21,1518) (AI(I),I=1,20),(KDSTL(I,2),I=1,20) 1 ,(RKDSTL(I,2),I=1,20) IF(NJ(2).EQ.1) GO TO 1800 DO 1781 J=NJ(2)+1,LJ-NJ(1) WRITE(21,1783) J-NJ(2) 1783 FORMAT(///5X,'POLAR ANGULAR DISTRIBUTION OF SPUTTERED PARTICLES ; 1LAYER 3 ; SPECIES ',I1//) WRITE(21,1518) (AI(I),I=1,20),(KDSTJ(I,J),I=1,20) 1 ,(RKDSTJ(I,J),I=1,20) 1781 CONTINUE 1800 CONTINUE c c hier wird der File for33 erzeugt c DO i=1,100 READ(33,'(A246)',ERR=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', 1 5x,'imp',2x,'backsc',3x,'trans',3x,'tried',4x,'negE',3x, 2 'impL1',3x,'impL2',3x,'impL3',3x,'impL4',3x,'impL5',3x,'impL6', 3 3x,'impL7',3x, 4 'range',6x,'straggeling',2x, 5 'Eback',7x,'sigEback',4x,'Etrans',6x,'SigEtrans',3x, 6 '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, 1 NH,IIM,IB,IT,tryE,negE, 2 (number_in_layer(k),k=1,7), 3 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 hier endet File for33 C C TOP AND FRONT LINES FOR MATRICES C C JE=DE C JA=DA C JG=DG C DO 32 J=2,NG1 C MAGB(J,1) = (J-1)*JG C MAGT(J,1) = (J-1)*JG C EMA(J,1)=DFLOTJ(J-1)*DG C EMAT(J,1)=DFLOTJ(J-1)*DG C 32 CONTINUE C DO 77 J=2,21 C MEAB(1,J) = J-1 C MEAT(1,J) = J-1 C MAGB(1,J) = J-1 C MAGT(1,J) = J-1 C EMA(1,J) = J-1 C EMAT(1,J) = J-1 C 77 CONTINUE C DO 1828 J=2,101 C MEAB(J,1) = J-1 C MEAT(J,1) = J-1 C MEPB(J,1) = J-1 C MEPB(1,J) = J-1 C MEPT(J,1) = J-1 C MEPT(1,J) = J-1 C 1828 CONTINUE C DO 1830 K=1,JT(3) C DO 1832 J=2,NG1 C MAGS(J,1,K) = (J-1)*JG C MAGST(J,1,K) = (J-1)*JG C MAGSA(J,1,K) = (J-1)*JG C 1832 CONTINUE C DO 1826 J=2,NA1 C MAGSA(1,J,K) = (J-1)*JA C 1826 CONTINUE C DO 1834 J=2,21 C MEAS(1,J,K) = J-1 C MEAST(1,J,K) = J-1 C MAGS(1,J,K) = J-1 C MAGST(1,J,K) = J-1 C 1834 CONTINUE C DO 1836 J=2,101 C MEAS(J,1,K) = J-1 C MEAST(J,1,K) = J-1 C 1836 CONTINUE C DO 1838 J=1,20 C MEASL(1,J,K)=J C MEASTL(1,J,K)=J C 1838 CONTINUE C DO 1841 IG2=1,NGIK,1 C DO 1843 J=2,21 C MEAGS(1,IG2,J,K) = J-1 C 1843 CONTINUE C DO 1845 J=2,101 C MEAGS(J,IG2,1,K) = J-1 C 1845 CONTINUE C 1841 CONTINUE C 1830 CONTINUE C DO 1840 IG2=1,NGIK,1 C DO 1842 J=2,21 C MEAGB(1,IG2,J) = J-1 C 1842 MEAGT(1,IG2,J) = J-1 C DO 1844 J=2,101 C MEAGB(J,IG2,1) = J-1 C 1844 MEAGT(J,IG2,1) = J-1 C 1840 CONTINUE C DO 1846 I=2,74 C 1846 ELOG(I)=10.D0**(I/12.D0)*10.D0**(-7.D0/6.D0) C TEMP=(10.D0**(1.D0/12.D0)-1.D0)*10.D0**(-7.D0/6.D0) C TEMPNH=TEMP*DFLOTJ(NH) C C BACKWARD SPUTTERING : MATRICES , ENERGY - ANGLE CORRELATIONS C C IF(ISPA.LT.10000) GO TO 1900 C DO 1850 J=1,JT(3) C EASL(2,J)=DFLOTJ(MEASL(2,21,J))/(DFLOTJ(NH)*0.1) C DO 1850 IESLOG=3,74 C 1850 EASL(IESLOG,J)=DFLOTJ(MEASL(IESLOG,21,J))/(TEMPNH* C 1 10.D0**((IESLOG-1)/12.D0)) C DO 1852 J=1,NJ(1) C WRITE(21,1854) J C 1854 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (BAC C 1KWARD SPUTTERED PARTICLES) ; 1. LAYER ; SPECIES',I2/) C do ima = 74,2,-1 C if(measl(ima,21,j).ne.0) goto 1855 C enddo C ima = 1 C 1855 ima = min(ima+2,74) C do ies = 1, ima C write (6, 1858) elog(ies), (measl(ies,ias,j),ias=1,21), C 1 easl(ies,j) C end do C write (6, 1858) elog(75), (measl(75,ias,j),ias=1,21), C 1 easl(75,j) cc DO 1856 IES=1,75 cc1856 WRITE(6,1858) ELOG(IES),(MEASL(IES,IAS,J),IAS=1,21),EASL(IES,J) C 1858 FORMAT(1X,1E12.4,20I5,I6,1E12.4) C WRITE(21,1884) J C 1884 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (BACKWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/) C do ima = 101,1,-1 C if(meas(ima,22,j).ne.0) goto 1883 C enddo C ima = 1 C 1883 ima = min(ima+2,101) C write (6, 1886) ((meas(iesp,iags,j),iags=1,22),iesp=1,ima) C write (6, 1886) (meas(102,iags,j),iags=1,22) C 1886 FORMAT(1X,I3,20I6,I8) Cc WRITE(6,1886) ((MEAS(IESP,IAGS,J),IAGS=1,22),IESP=1,102) Cc1886 FORMAT(1X,21I5,I7) C IF(ALPHA.LT.1.) GO TO 1878 C DO 1870 IG2=1,NGIK,1 C EEE = IG2*DGI C WRITE(21,1872) EEE,J C 1872 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) AT AZIMUTHAL ANGLE =',F5.1,' (BACKWARD SPUTTERED ATOMS) , 1.LAYE C 2R , SPECIES',I2/) C do ima = 101,1,-1 C if(meags(ima,ig2,22,j).ne.0) goto 1885 C enddo C ima = 1 C 1885 ima = min(ima+2,101) C do iesp = 1, ima C write (6, 1886) (meags(iesp,ig2,iags,j),iags=1,22) C end do C write (6, 1886) (meags(102,ig2,iags,j),iags=1,22) Cc DO 1870 IE=1,102 Cc WRITE(6,1886) (MEAGS(IE,IG2,IAGS,J),IAGS=1,22) C 1870 CONTINUE C WRITE(21,1889) J C 1889 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN DEGREES (BACKWARD C 1 SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/) C WRITE(21,1887) ((MAGSA(IG,IA,J),IA=1,32),IG=1,62) C 1887 FORMAT(1X,31I4,I6) C 1878 CONTINUE C WRITE(21,1888) J C 1888 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (BACKWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/) C WRITE(21,1886) ((MAGS(IG,IAGS,J),IAGS=1,22),IG=1,62) C 1852 CONTINUE C IF(L.EQ.1) GO TO 1900 C if(ispal(2).eq.0) goto 1900 C DO 1862 J=NJ(1)+1,JT(3) C WRITE(21,1864) J-NJ(1) C 1864 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (BAC C 1KWARD SPUTTERED PARTICLES) , 2. LAYER , SPECIES',I2/) C do ima = 74,1,-1 C if(measl(ima,21,j).ne.0) goto 1865 C enddo C ima = 1 C 1865 ima = min(ima+2,74) C do ies = 1, ima C write (6, 1858) elog(ies), (measl(ies,ias,j),ias=1,21) C 1 , easl(ies,j) C end do C write (6, 1858) elog(75), (measl(75,ias,j),ias=1,21) C 1 , easl(75,j) Cc DO 1866 IES=1,75 Cc1866 WRITE(6,1858) ELOG(IES),(MEASL(IES,IAS,J),IAS=1,21),EASL(IES,J) C WRITE(21,1894) J-NJ(1) C 1894 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (BACKWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/) C do ima = 101,1,-1 C if(meas(ima,22,j).ne.0) goto 1895 C enddo C ima = 1 C 1895 ima = min(ima+2,101) C WRITE(21,1886)((meas(iesp,iags,j),iags=1,22),iesp=1,ima) C WRITE(21,1886)(meas(102,iags,j),iags=1,22) Cc WRITE(6,1886) ((MEAS(IESP,IAGS,J),IAGS=1,22),IESP=1,102) C WRITE(21,1898) J-NJ(1) C 1898 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (BACKWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/) C WRITE(21,1886) ((MAGS(IG,IAGS,J),IAGS=1,22),IG=1,62) C 1862 CONTINUE C 1900 CONTINUE CC CC FORWARD SPUTTERING : MATRICES , ENERGY - ANGLE CORRELATIONS CC C IF(ISPAT.LT.10000) GO TO 2000 C JTJ=JT(3) C IF(L.EQ.3) JTJ=LJ-NJ(1) C DO 1950 J=1,JTJ C EASTL(2,J)=DFLOTJ(MEASTL(2,21,J))/(DFLOTJ(NH)*0.1D0) C DO 1950 IESLOG=3,74 C 1950 EASTL(IESLOG,J)=DFLOTJ(MEASTL(IESLOG,21,J))/(TEMPNH* C 1 10.D0**((IESLOG-1)/12.D0)) C IF(L.EQ.3) GO TO 1953 C DO 1952 J=1,NJ(1) C WRITE(21,1954) J C 1954 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (FOR C 1WARD SPUTTERED PARTICLES) , 1. LAYER , SPECIES',I2/) C do ima = 74,2,-1 C if(meastl(ima,21,j).ne.0) goto 1955 C enddo C ima = 1 C 1955 ima = min(ima+2,74) C do ies = 1, ima C write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21), C 1 eastl(ies,j) C end do C write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21), C 1 eastl(75,j) Cc DO 1956 IES=1,75 Cc1956 WRITE(6,1858) ELOG(IES),(MEASTL(IES,IAS,J),IAS=1,21),EASTL(IES,J) C WRITE(21,1984) J C 1984 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (FORWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/) C do ima = 101,1,-1 C if(meast(ima,22,j).ne.0) goto 1983 C enddo C ima = 1 C 1983 ima = min(ima+2,101) C write (6, 1886) ((meast(iesp,iags,j),iags=1,22),iesp=1,ima) C write (6, 1886) (meast(102,iags,j),iags=1,22) Cc WRITE(6,1886) ((MEAST(IESP,IAGS,J),IAGS=1,22),IESP=1,102) C WRITE(21,1988) J C 1988 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (FORWARD SPUTTERED PARTICLES) , 1.LAYER , SPECIES',I2/) C WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62) C 1952 CONTINUE C 1953 CONTINUE C IF(L.EQ.1) GO TO 2000 C IF(L.EQ.3) GO TO 1961 C JTK=NJ(1)+1 C JTL=JT(3) C GO TO 1963 C 1961 JTK=1 C JTL=NJ(2) C 1963 DO 1962 J=JTK,JTL C WRITE(21,1964) J-JTK+1 C 1964 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (FOR C 1WARD SPUTTERED PARTICLES) , 2. LAYER , SPECIES',I2/) C do ima = 74,1,-1 C if(meastl(ima,21,j).ne.0) goto 1965 C enddo C ima = 1 C 1965 ima = min(ima+2,74) C do ies = 1, ima C write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21) C 1 , eastl(ies,j) C end do C write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21) C 1 , eastl(75,j) Cc DO 1966 IES=1,75 Cc1966 WRITE(6,1858) ELOG(IES),(MEASTL(IES,IAS,J),IAS=1,21),EASTL(IES,J) C WRITE(21,1994) J-JTK+1 C 1994 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (FORWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/) C do ima = 101,1,-1 C if(meast(ima,22,j).ne.0) goto 1995 C enddo C ima = 1 C 1995 ima = min(ima+2,101) C WRITE(21,1886)((meast(iesp,iags,j),iags=1,22),iesp=1,ima) C WRITE(21,1886)(meast(102,iags,j),iags=1,22) Cc WRITE(6,1886) ((MEAST(IESP,IAGS,J),IAGS=1,22),IESP=1,102) C WRITE(21,1998) J-JTK+1 C 1998 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (FORWARD SPUTTERED PARTICLES) , 2.LAYER , SPECIES',I2/) C WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62) C 1962 CONTINUE C IF(L.LT.3) GO TO 2000 C DO 1972 J=NJ(2)+1,LJ-NJ(1) C WRITE(21,1974) J-NJ(2) C 1974 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (FOR C 1WARD SPUTTERED PARTICLES) , 3. LAYER , SPECIES',I2/) C do ima = 74,1,-1 C if(meastl(ima,21,j).ne.0) goto 1973 C enddo C ima = 1 C 1973 ima = min(ima+2,74) C do ies = 1, ima C write (6, 1858) elog(ies), (meastl(ies,ias,j),ias=1,21) C 1 , eastl(ies,j) C end do C write (6, 1858) elog(75), (meastl(75,ias,j),ias=1,21) C 1 , eastl(75,j) Cc DO 1976 IES=1,75 Cc1976 WRITE(6,1858) ELOG(IES),(MEASTL(IES,IAS,J),IAS=1,21),EASTL(IES,J) C WRITE(21,1975) J-NJ(2) C 1975 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (FORWARD SPUTTERED PARTICLES) , 3.LAYER , SPECIES',I2/) C do ima = 101,1,-1 C if(meast(ima,22,j).ne.0) goto 1977 C enddo C ima = 1 C 1977 ima = min(ima+2,101) C WRITE(21,1886)((meast(iesp,iags,j),iags=1,22),iesp=1,ima) C WRITE(21,1886)(meast(102,iags,j),iags=1,22) Cc WRITE(6,1886) ((MEAST(IESP,IAGS,J),IAGS=1,22),IESP=1,102) C WRITE(21,1978) J-NJ(2) C 1978 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (FORWARD SPUTTERED PARTICLES) , 3.LAYER , SPECIES',I2/) C WRITE(21,1886) ((MAGST(IG,IAGS,J),IAGS=1,22),IG=1,62) C 1972 CONTINUE CC DO 34 IG2=1,NGIK,1 CC EEE = IG2*DGI CC WRITE(6,912) EEE CC 912 FORMAT(1H1,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.0 CC 15) AT AZIMUTHAL ANGLE =',F5.1,' (SPUTTERED ATOMS)'//) CC DO 42 IE=2,101 CC 42 MEAGS(102,IG2,22) = MEAGS(102,IG2,22)+MEAGS(IE,IG2,22) CC DO 34 IE=1,102 CC WRITE(6,980) (MEAGS(IE,IG2,IAGS),IAGS=1,22) CC 34 CONTINUE CC IF(ALPHA.LT.1.) GO TO 8009 CC DO 8001 IG3=1,NGIK,1 CC EE1 = IG3*DGI CC WRITE(6,8002) EE1 CC8002 FORMAT(1H1,' LOG ENERGY - POLAR ANGLE IN COS-INTERVALS (0.05) AT CC 1 AZIMUTHAL ANGLE =',F5.1,' (SPUTTERED ATOMS)'//) CC DO 8003 J=1,20 CC8003 MEAGSL(1,IG3,J)=J CC IF(MEAGS(102,IG3,22).EQ.0) MEAGS(102,IG3,22)=1 CC EAGSL(2)=DFLOAT(MEAGSL(2,IG3,21))/(DFLOAT(MEAGS(102,IG3,22))*0.1) CC DO 8004 IESLOG=3,74 CC8004 EAGSL(IESLOG)=DFLOAT(MEAGSL(IESLOG,IG3,21))/(DFLOAT(MEAGS(102,IG3,22 CC ?))*TEMP*10.**((IESLOG-1)/12.)) CC DO 8005 IES=1,75 CC8005 WRITE(6,8600) ELOG(IES),(MEAGSL(IES,IG3,IAS),IAS=1,21),EAGSL(IES) CC8001 CONTINUE C 2000 CONTINUE CC CC BACKSCATTERING : MATRICES , ENERGY - ANGULAR CORRELATIONS CC C IF(IB.LT.10000) GO TO 2100 C DO 2002 J=1,20 C 2002 MEABL(1,J)=J C EABL(2)=DFLOTJ(MEABL(2,21))/(DFLOTJ(NH)*0.1D0) C DO 2004 IERLOG=3,74 C 2004 EABL(IERLOG)=DFLOTJ(MEABL(IERLOG,21))/(TEMPNH* C #10.D0**((IERLOG-1)/12.D0)) C WRITE(21,2006) C 2006 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (BAC C 1KSCATTERED PROJECTILES)'/) C do ima = 74,1,-1 C if(meabl(ima,21).ne.0) goto 2005 C enddo C ima = 1 C 2005 ima = min(ima+2,74) C do ies = 1, ima C WRITE(21,1858)elog(ies),(meabl(ies,iag),iag=1,21),eabl(ies) C end do C WRITE(21,1858)elog(75),(meabl(75,iag),iag=1,21),eabl(75) Cc DO 2008 IES=1,75 Cc2008 WRITE(6,1858) ELOG(IES),(MEABL(IES,IAG),IAG=1,21),EABL(IES) C IF(ALPHA.LT.1.) GO TO 2010 C DO 2012 IG2=1,NGIK,1 C EEE = IG2*DGI C WRITE(21,2014) EEE C 2014 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) AT AZIMUTHAL ANGLE =',F5.1,' (BACKSCATTERED PROJECTILES)'/) C do ima = 101,1,-1 C if(meagb(ima,ig2,22).ne.0) goto 2015 C enddo C ima = 1 C 2015 ima = min(ima+2,101) C write (6, 1886) ((meagb(ie,ig2,iagb),iagb=1,22),ie=1,ima) C write (6, 1886) (meagb(102,ig2,iagb),iagb=1,22) Cc2012 WRITE(6,1886) ((MEAGB(IE,IG2,IAGB),IAGB=1,22),IE=1,102) C 2012 continue C 2010 CONTINUE C IF(E0.LT.0.) GO TO 2052 C WRITE(21,2016) C 2016 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (BACKSCATTERED PROJECTILES)'/) C GO TO 2054 C 2052 WRITE(21,2056) C 2056 FORMAT(//,' ENERGY(E IN 0.1*TI) - POLAR ANGLE IN COS-INTERVALS (0. C 105) (BACKSCATTERED PROJECTILES)'/) C do ima = 101,1,-1 C if(meab(ima,22).ne.0) goto 2017 C enddo C ima = 1 C 2017 ima = min(ima+2,101) C write (6, 1886) ((meab(ie,iagb),iagb=1,22),ie=1,ima) C write (6, 1886) (meab(102,iagb),iagb=1,22) Cc2054 WRITE(6,1886) ((MEAB(IE,IAGB),IAGB=1,22),IE=1,NE) C 2054 continue C WRITE(21,2018) C 2018 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (BACKSCATTERED PROJECTILES)'/) C WRITE(21,1886) ((MAGB(IG,IAGB),IAGB=1,22),IG=1,62) C WRITE(21,2022) C 2022 FORMAT(//,' AZIMUTH.ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (BACKSCATTERED ENERGY)'/) C WRITE(21,2027) (EMA(01,IAGB),IAGB=1,11) C WRITE(21,2025) ((EMA(IG,IAGB),IAGB=1,11),IG=2,NG) C WRITE(21,2028) C WRITE(21,2031) EMA(1,1),(EMA(1,IAGB),IAGB=12,22) C DO 2029 IG=2,NG C WRITE(21,2026) EMA(IG,1),(EMA(IG,IAGB),IAGB=12,22) C 2029 CONTINUE C 2025 FORMAT(1X,1F5.0,10E11.4) C 2026 FORMAT(1X,1F5.0,11E11.4) C 2027 FORMAT(1X,1F5.0,10F11.0) C 2028 FORMAT(/) C 2031 FORMAT(1H1,1X,1F5.0,11F11.0) CC IF(E0.LT.0.) GO TO 2058 CC WRITE(6,2032) CC2032 FORMAT(1H1,1X,'ENERGY(IN % OF E0) - PATHLENGTH(IN UNITS OF CW) CC 1 (BACKSCATTERED PROJECTILES)'/) CC GO TO 2060 CC2058 WRITE(6,2062) CC2062 FORMAT(1H1,1X,'ENERGY(E IN 0.1*TI) - PATHLENGTH(IN UNITS OF CW) CC 1 (BACKSCATTERED PROJECTILES)'/) CC2060 DO 2034 II=1,3 CC INE=II*25+1 CC INA=INE-24 CC DO 2040 IE=1,NE CC WRITE(6,2036) MEPB(IE,1),(MEPB(IE,IPB),IPB=INA,INE) CC2040 CONTINUE CC WRITE(6,2028) CC2034 CONTINUE CC DO 2042 IE=1,NE CC WRITE(6,2038) MEPB(IE,1),(MEPB(IE,IPB),IPB=77,102) CC2042 CONTINUE C 2036 FORMAT(1X,26I4) C 2038 FORMAT(1X,26I4,I6) C 2100 CONTINUE CC CC TRANSMISSION : MATRICES , ENERGY - ANGULAR CORRELATIONS CC C IF(IT.LT.10000) GO TO 9000 C DO 2102 J=1,20 C 2102 MEATL(1,J)=J C EATL(2)=DFLOTJ(MEATL(2,21))/(DFLOTJ(NH)*0.1D0) C DO 2104 IERLOG=3,74 C 2104 EATL(IERLOG)=DFLOTJ(MEATL(IERLOG,21))/(TEMPNH* C 1 10.D0**((IERLOG-1)/12.D0)) C WRITE(21,2106) C 2106 FORMAT(//,' LOG ENERGY - COS OF EMISSION ANGLE (0.05 STEPS) (TRA C 1NSMITTED PROJECTILES)'/) C do ima = 74,1,-1 C if(meatl(ima,21).ne.0) goto 2105 C enddo C ima = 1 C 2105 ima = min(ima+2,74) C do ies = 1, ima C WRITE(21,1858)elog(ies),(meatl(ies,iag),iag=1,21),eatl(ies) C end do C WRITE(21,1858)elog(75),(meatl(75,iag),iag=1,21),eatl(75) Cc DO 2108 IES=1,75 Cc2108 WRITE(21,1858) ELOG(IES),(MEATL(IES,IAG),IAG=1,21),EATL(IES) C IF(ALPHA.LT.1.) GO TO 2110 C DO 2112 IG2=1,NGIK,1 C EEE = IG2*DGI C WRITE(21,2114) EEE C 2114 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) AT AZIMUTHAL ANGLE =',F5.1,' (TRANSMITTED PROJECTILES)'/) C do ima = 101,1,-1 C if(meagt(ima,ig2,22).ne.0) goto 2115 C enddo C ima = 1 C 2115 ima = min(ima+2,101) C write (21,1886) ((meagt(ie,ig2,iagb),iagb=1,22),ie=1,ima) C write (21,1886) (meagt(102,ig2,iagb),iagb=1,22) Cc2112 WRITE(6,1886) ((MEAGT(IE,IG2,IAGB),IAGB=1,22),IE=1,102) C 2112 continue C 2110 CONTINUE C WRITE(21,2116) C 2116 FORMAT(//,' ENERGY(E/E0 IN %) - POLAR ANGLE IN COS-INTERVALS (0.05 C 1) (TRANSMITTED PROJECTILES)'/) C do ima = 101,1,-1 C if(meat(ima,22).ne.0) goto 2117 C enddo C ima = 1 C 2117 ima = min(ima+2,101) C write (6, 1886) ((meat(ie,iagb),iagb=1,22),ie=1,ima) C write (6, 1886) (meat(102,iagb),iagb=1,22) Cc WRITE(6,1886) ((MEAT(IE,IAGB),IAGB=1,22),IE=1,NE) C WRITE(21,2118) C 2118 FORMAT(//,' AZIMUTHAL ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (TRANSMITTED PROJECTILES)'/) C WRITE(21,1886) ((MAGT(IG,IAGB),IAGB=1,22),IG=1,62) C WRITE(21,2122) C 2122 FORMAT(//,' AZIMUTH.ANGLE - POLAR ANGLE IN COS-INTERVALS (0.05) C 1 (TRANSMITTED ENERGY)'/) C WRITE(21,2127) (EMAT(01,IAGB),IAGB=1,11) C WRITE(21,2125) ((EMAT(IG,IAGB),IAGB=1,11),IG=2,NG) C WRITE(21,2028) C WRITE(21,2131) EMAT(1,1),(EMAT(1,IAGB),IAGB=12,22) C DO 2129 IG=2,NG C WRITE(21,2126) EMAT(IG,1),(EMAT(IG,IAGB),IAGB=12,22) C 2129 CONTINUE C 2125 FORMAT(1X,1F5.0,10E11.4) C 2126 FORMAT(1X,1F5.0,11E11.4) C 2127 FORMAT(1X,1F5.0,10F11.0) C 2131 FORMAT(1H1,1X,1F5.0,11F11.0) C GO TO 9000 CC WRITE(6,2132) CC2132 FORMAT(1H1,1X,'ENERGY(IN % OF E0) - PATHLENGTH(IN UNITS OF CW) CC 1 (TRANSMITTED PROJECTILES)'/) CC DO 2134 II=1,3 CC INE=II*25+1 CC INA=INE-24 CC DO 2140 IE=1,NE CC WRITE(6,2036) MEPT(IE,1),(MEPT(IE,IPT),IPT=INA,INE) CC2140 CONTINUE CC WRITE(6,2028) CC2134 CONTINUE CC DO 2142 IE=1,NE CC WRITE(6,2038) MEPT(IE,1),(MEPT(IE,IPT),IPT=77,102) CC2142 CONTINUE C 9000 CONTINUE CC WRITE(6,*) INEL,L,LJ CC DO 9875 J=1,180 CC IANF=J*7-6 CC IEND=(J+1)*7-7 CC WRITE(6,9876) (ESVDL(I),I=IANF,IEND) CC9876 FORMAT(1X,7E11.4) CC9875 CONTINUE CC CC DATA ON DISC CC c WRITE(17) Z1,M1,E0,ALPHA,EF,ESB,SHEATH c 1 ,NH,RI,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 c WRITE(17) (DX(I),I=1,3),(RHO(I),I=1,3),(CK(I),I=1,3) c 1 ,((ZT(I,J),J=1,5),I=1,3),((MT(I,J),J=1,5),I=1,3) c 2 ,((CO(I,J),J=1,5),I=1,3),((SBE(I,J),J=1,5),I=1,3) c 3 ,((ED(I,J),J=1,5),I=1,3),((BE(I,J),J=1,5),I=1,3) c WRITE(17) TI,ZARG,VELC c 1 ,HLM,HLMT,SU,SUT,XC,RT,INEL,L,LJ c 2 ,NPROJ,KIB,KIT,MAXA,NALL,NPA,NSA,KIS,KIST c 3 ,IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,ISPAT,ESPAT c 4 ,FIX0,SEX,THX,FOX,SIGMAX,DFIX0,DSEX,DTHX c 5 ,FIR0,SER,THR,FOR,SIGMAR,DFIR0,DSER,DTHR c 6 ,FIP0,SEP,THP,FOP,SIGMAP,DFIP0,DSEP,DTHP c 7 ,AVNLI,VANLI,SIGNLI,DFINLI c 8 ,AVILI,VAILI,SIGILI,DFIILI c WRITE(17) AVCSUM,AVCDIS c 1 ,FIE0,SEE,THE,FOE,SIGMAE,DFIE0,DSEE,DTHE c 2 ,FIW0,SEW,THW,FOW,SIGMAW,DFIW0,DSEW,DTHW c 3 ,FII0,SEI,THI,FOI,SIGMAI,DFII0,DSEI,DTHI c 4 ,FIS0,SES,THS,FOS,SIGMAS,DFIS0,DSES,DTHS c 5 ,IIRP,TRIRP,IIPL,TION,TDMGN,TCASMO,TPHON,TDENT c WRITE(17) RN,RE,EMEANR,EMEAN,TN,TE,TMEANR,EMEANT c 1 ,FIB0,SEB,THB,FOB,SIGMAB,DFIB0,DSEB,DTHB c 2 ,FIPB0,SEPB,THPB,FOPB,SIGMPB,DFIPB0,DSEPB,DTHPB c 3 ,AVNLB,VANLB,SIGNLB,DFINLB c 4 ,AVILB,VAILB,SIGILB,DFIILB c WRITE(17) FIT0,SET,THT,FOT,SIGMAT,DFIT0,DSET,DTHT c 1 ,FIPT0,SEPT,THPT,FOPT,SIGMPT,DFIPT0,DSEPT,DTHPT c 2 ,AVNLT,VANLT,SIGNLT,DFINLT c 3 ,AVILT,VAILT,SIGILT,DFIILT c WRITE(17) (IRP(I),I=0,100),(RIRP(I),I=0,100) c 1 ,(IPL(I),I=1,100),(ION(I),I=1,100),(DMGN(I),I=1,100) c 2 ,(CASMOT(I),I=1,100),(PHON(I),I=1,100),(DENT(I),I=1,100) c WRITE(17) (FIESB(J),J=1,10),(SEESB(J),J=1,10),(THESB(J),J=1,10) c 1 ,(FOESB(J),J=1,10),(SGMESB(J),J=1,10) c 2 ,(DFIESB(J),J=1,10),(DSEESB(J),J=1,10) c 3 ,(DTHESB(J),J=1,10) c WRITE(17) ((ELE(I,J),J=1,15),I=1,100),((ELI(I,J),J=1,15),I=1,100) c 1 ,((ELP(I,J),J=1,15),I=1,100) c 2 ,(ELET(J),J=1,15),(ELIT(J),J=1,15),(ELPT(J),J=1,15) c WRITE(17) (AI(I),I=1,20),(KADB(I),I=1,20),(KADT(I),I=1,20) c 1 ,(RKADB(I),I=1,20),(RKADT(I),I=1,20) c WRITE(17) (KADS(I),I=1,20),(KADST(I),I=1,20) c 1 ,(RKADS(I),I=1,20),(RKADST(I),I=1,20) c WRITE(17) ((KADRIP(I,J),J=1,10),I=1,20) c 1 ,((KADRIS(I,J),J=1,10),I=1,20) c 2 ,((KADROP(I,J),J=1,10),I=1,20) c 3 ,((KADROS(I,J),J=1,10),I=1,20) cC 4 ,RKDRIP(20),RKDRIS(20),RKDROP(20),RKDROS(20) c WRITE(17) ((KADSJ(I,J),J=1,10),I=1,20) c 1 ,((RKADSJ(I,J),J=1,10),I=1,20) c 2 ,((KADSL(I,J),J=1,2),I=1,20) c 3 ,((RKADSL(I,J),J=1,2),I=1,20) c WRITE(17) ((KDSTJ(I,J),J=1,10),I=1,20) c 1 ,((RKDSTJ(I,J),J=1,10),I=1,20) c 2 ,((KDSTL(I,J),J=1,2),I=1,20) c 3 ,((RKDSTL(I,J),J=1,2),I=1,20) c WRITE(17) (IBSP(I),I=1,15),(EBSP(I),I=1,15) c 1 ,(SPY(I),I=1,15),(SPE(I),I=1,15) c 2 ,(REY(I),I=1,15),(EMSP(I),I=1,15) c 3 ,(ISPAL(I),I=1,3),(ESPAL(I),I=1,3) c WRITE(17) (ISPIP(I),I=1,15),(ISPIS(I),I=1,15) c 1 ,(ISPOP(I),I=1,15),(ISPOS(I),I=1,15) c 2 ,(ESPIP(I),I=1,15),(ESPIS(I),I=1,15) c 3 ,(ESPOP(I),I=1,15),(ESPOS(I),I=1,15) c 4 ,(RIP(I),I=1,15),(RIS(I),I=1,15) c 5 ,(ROP(I),I=1,15),(ROS(I),I=1,15) c 6 ,(REIP(I),I=1,15),(REIS(I),I=1,15) c 7 ,(REOP(I),I=1,15),(REOS(I),I=1,15) c WRITE(17) (ITSP(I),I=1,15),(ETSP(I),I=1,15) c 1 ,(SPYT(I),I=1,15),(SPET(I),I=1,15) c 2 ,(REYT(I),I=1,15),(EMSPT(I),I=1,15) c 3 ,(ISPALT(I),I=1,3),(ESPALT(I),I=1,3) c WRITE(17) (ISPIPT(I),I=1,15),(ISPIST(I),I=1,15) c 1 ,(ISPOPT(I),I=1,15),(ISPOST(I),I=1,15) c 2 ,(ESPIPT(I),I=1,15),(ESPIST(I),I=1,15) c 3 ,(ESPOPT(I),I=1,15),(ESPOST(I),I=1,15) c 4 ,(RIPT(I),I=1,15),(RIST(I),I=1,15) c 5 ,(ROPT(I),I=1,15),(ROST(I),I=1,15) c 6 ,(REIPT(I),I=1,15),(REIST(I),I=1,15) c 7 ,(REOPT(I),I=1,15),(REOST(I),I=1,15) c WRITE(17) ((MEAB(I,J),J=1,22),I=1,102) c 1 ,((MAGB(I,J),J=1,22),I=1,62) c 2 ,(((MEAGB(I,J,K),K=1,22),J=1,36),I=1,102) c 3 ,((EMA(I,J),J=1,22),I=1,62),(ELOG(I),I=1,75) c 4 ,(EABL(I),I=1,75),((MEABL(I,J),J=1,21),I=1,75) c 5 ,((MEPB(I,J),J=1,102),I=1,102) c WRITE(17) ((MEAT(I,J),J=1,22),I=1,102) c 1 ,((MAGT(I,J),J=1,22),I=1,62) c 2 ,(((MEAGT(I,J,K),K=1,22),J=1,36),I=1,102) c 3 ,((EMAT(I,J),J=1,22),I=1,62) c 4 ,(EATL(I),I=1,75),((MEATL(I,J),J=1,21),I=1,75) c 5 ,((MEPT(I,J),J=1,102),I=1,102) c WRITE(17) (((MEAS(I,J,K),K=1,10),J=1,22),I=1,102) c 1 ,(((MAGS(I,J,K),K=1,10),J=1,22),I=1,62) c 2 ,((EASL(I,J),J=1,10),I=1,75) c 3 ,(((MEASL(I,J,K),K=1,10),J=1,21),I=1,75) c WRITE(17) (((MEAST(I,J,K),K=1,10),J=1,22),I=1,102) c 1 ,(((MAGST(I,J,K),K=1,10),J=1,22),I=1,62) c 2 ,((EASTL(I,J),J=1,10),I=1,75) c 3 ,(((MEASTL(I,J,K),K=1,10),J=1,21),I=1,75) c WRITE(17) ((((MEAGS(I,J,K,MN),MN=1,10),K=1,22),J=1,12),I=1,102) c 1 ,(((MAGSA(I,J,K),K=1,10),J=1,32),I=1,62) CC 1 ,((((MEAGST(I,J,K,L),L=1,10),K=1,22),J=1,36),I=1,102) c WRITE(17) ((ELD(I,J),I=1,100),J=1,15) c WRITE(17) XSUM,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM c WRITE(17) EB,EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM c 1 ,EB1SUL,EB2SUL,EB3SUL,EB4SUL,EB5SUL,EB6SUL c WRITE(17) (EBSP(J),J=1,15),(SPE2S(J),J=1,15),(SPE3S(J),J=1,15) c 1 ,(SPE4S(J),J=1,15),(SPE5S(J),J=1,15),(SPE6S(J),J=1,15) c WRITE(17) (SPE1SL(J),J=1,15),(SPE2SL(J),J=1,15),(SPE3SL(J),J=1,15) c 1 ,(SPE4SL(J),J=1,15),(SPE5SL(J),J=1,15) c 2 ,(SPE6SL(J),J=1,15) c WRITE(17) ((ICD(I,J),J=1,15),I=1,100),((ICDR(I,J),J=1,15),I=1,100) c WRITE(17) (((ICDIRI(I,J,K),K=1,15),J=1,15),I=1,100) c 1 ,((ICDIRN(I,J),J=1,15),I=1,100) c write(17) exi1s,exi2s,exi3s,exi4s,exi5s,exi6s c 1 ,coss1s,coss2s,coss3s,coss4s,coss5s,coss6s c write(17) ibl,(ibsp(i),i=1,15) C CLOSE(UNIT=21) CLOSE(UNIT=22) CLOSE(UNIT=99) 8000 STOP END C C SUBROUTINE MAGICKRC(C2,S2,B,R,EPS,N) C DIMENSION C2(N),S2(N),B(N),R(N),EPS(N),V(N),V1(N),TEST(N) C DIMENSION EX1(N),EX2(N),EX3(N) C IVMIN=1 C IVMAX=N C C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.) C C DO 105 IV=IVMIN,IVMAX C KRYPTON-CARBON POTENTIAL C EX1(IV)=DEXP(-.278544*R(IV)) C EX2(IV)=DEXP(-.637174*R(IV)) C EX3(IV)=DEXP(-1.919249*R(IV)) C RR1=1./R(IV) C V(IV)=(.190945*EX1(IV)+.473674*EX2(IV)+.335381*EX3(IV))*RR1 C V1(IV)=-(V(IV)+.053186584080*EX1(IV)+.301812757276*EX2(IV)+ C 1 .643679648869*EX3(IV))*RR1 C FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) C FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1. C Q=FR/FR1 C R(IV)=R(IV)-Q C TEST(IV)=DABS(Q/R(IV)).GT.0.001 C 105 CONTINUE C GET MAX AND MIN INDEX OF TEST FAILURES C IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) C IF(IVMIN.GT.IVMAX) GO TO 106 C IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) C IF(IVMIN.GT.IVMAX) GO TO 106 C GO TO 104 C 106 DO 108 IV=1,IH1 C ROCINV=-0.5*V1(IV)/(EPS(IV)-V(IV)) C SQE=DSQRT(EPS(IV)) C CC=(.235800+SQE)/(.126000+SQE) C AA=2.*EPS(IV)*(1.+(1.0144/SQE))*B(IV)**CC C FF=(DSQRT(AA*AA+1.)-AA)*((69350.+EPS(IV))/(83550.+EPS(IV))) C DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.) C C=(ROCINV*(B(IV)+DELTA)+1.)/(ROCINV*R(IV)+1.) C C2(IV)=DMIN1(1.0,C*C) C 108 S2(IV)=1.-C2(IV) C RETURN C END C C SUBROUTINE MAGICMOL(C2,S2,B,R,EPS,N) C DIMENSION C2(N),S2(N),B(N),R(N),EPS(N),V(N),V1(N),TEST(N) C DIMENSION EX1(N),EX2(N),EX3(N) C IVMIN=1 C IVMAX=N C C MAGIC (DETERMINATION OF SCATTERING ANGLE : MOLIERE POT.) C C DO 105 IV=IVMIN,IVMAX C MOLIERE POTENTIAL C EX1(IV)=DEXP(-.3*R(IV)) C EX2(IV)=DEXP(-1.2*R(IV)) C EX3(IV)=DEXP(-6.0*R(IV)) C RR1=1./R(IV) C V(IV)=(.35*EX1(IV)+.55*EX2(IV)+.10*EX3(IV))*RR1 C V1(IV)=-(V(IV)+.105*EX1(IV)+.66*EX2(IV)+.6*EX3(IV))*RR1 C FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) C FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1. C Q=FR/FR1 C R(IV)=R(IV)-Q C TEST(IV)=DABS(Q/R(IV)).GT.0.001 C 105 CONTINUE C GET MAX AND MIN INDEX OF TEST FAILURES C IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) C IF(IVMIN.GT.IVMAX) GO TO 106 C IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) C IF(IVMIN.GT.IVMAX) GO TO 106 C GO TO 104 C 106 DO 108 IV=1,IH1 C ROCINV=-0.5*V1(IV)/(EPS(IV)-V(IV)) C SQE=DSQRT(EPS(IV)) C CC=(.009611+SQE)/(.005175+SQE) C AA=2.*EPS(IV)*(1.+(0.6743/SQE))*B(IV)**CC C FF=(DSQRT(AA*AA+1.)-AA)*((6.314+EPS(IV))/(10.+EPS(IV))) C DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.) C C=(ROCINV*(B(IV)+DELTA)+1.)/(ROCINV*R(IV)+1.) C C2(IV)=DMIN1(1.0,C*C) C 108 S2(IV)=1.-C2(IV) C RETURN C END C C SUBROUTINE MAGICZBL(C2,S2,B,R,EPS,N) C DIMENSION C2(N),S2(N),B(N),R(N),EPS(N),V(N),V1(N),TEST(N) C DIMENSION EX1(N),EX2(N),EX3(N),EX4(N) C IVMIN=1 C IVMAX=N C C MAGIC (DETERMINATION OF SCATTERING ANGLE : ZBL POT.) C C DO 105 IV=IVMIN,IVMAX C ZBL POTENTIAL C EX1(IV)=DEXP(-.20162*R(IV)) C EX2(IV)=DEXP(-.4029*R(IV)) C EX3(IV)=DEXP(-.94229*R(IV)) C EX4(IV)=DEXP(-3.1998*R(IV)) C RR1=1./R(IV) C V(IV)=(.02817*EX1(IV)+.28022*EX2(IV)+.50986*EX3(IV)+ C 1 .18175*EX4(IV))*RR1 C V1(IV)=-(V(IV)+.0056796354*EX1(IV)+.112900638*EX2(IV)+ C 1 .4804359794*EX3(IV)+.581563650*EX4(IV))*RR1 C FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) C FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1. C Q=FR/FR1 C R(IV)=R(IV)-Q C TEST(IV)=DABS(Q/R(IV)).GT.0.001 C 105 CONTINUE C GET MAX AND MIN INDEX OF TEST FAILURES C IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) C IF(IVMIN.GT.IVMAX) GO TO 106 C IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) C IF(IVMIN.GT.IVMAX) GO TO 106 C GO TO 104 C 106 DO 108 IV=1,IH1 C ROCINV=-0.5*V1(IV)/(EPS(IV)-V(IV)) C SQE=DSQRT(EPS(IV)) C CC=(.011615+SQE)/(.0071222+SQE) C AA=2.*EPS(IV)*(1.+(0.99229/SQE))*B(IV)**CC C FF=(DSQRT(AA*AA+1.)-AA)*((9.3066+EPS(IV))/(14.813+EPS(IV))) C DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.) C C=(ROCINV*(B(IV)+DELTA)+1.)/(ROCINV*R(IV)+1.) C C2(IV)=DMIN1(1.0,C*C) C 108 S2(IV)=1.-C2(IV) C RETURN C END C SUBROUTINE MOMENTS(FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM, # X1S,X2S,X3S,X4S,X5S,X6S,Y) !DEC$REAL:8 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 C C IF(Y.EQ.0.D0.OR.Y.EQ.1.D0) GO TO 10 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 C 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) !DEC$REAL:8 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 C IF(Y.EQ.0.D0.OR.Y.EQ.1.D0) GO TO 10 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 C SUBROUTINE MOMENT(X1SY,X2SY,X3SY,X4SY,X5SY,X6SY # ,X1S,X2S,X3S,X4S,X5S,X6S,Y) !DEC$REAL:8 REAL*8 X1SY,X2SY,X3SY,X4SY,X5SY,X6SY,X1S,X2S,X3S,X4S,X5S,X6S,Y LOGICAL EQUAL C IF(Y.EQ.0.0D0) GO TO 10 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 C SUBROUTINE DIRCOS(COSX,COSY,COSZ,SINE,CPSI,SPSI,CPHI,SPHI,N) !DEC$REAL:8 INTEGER N,IV REAL*8 COSX(N),COSY(N),COSZ(N),SINE(N),CPSI(N),SPSI(N) # ,CPHI(N),SPHI(N) REAL*8 SRAT,CX2,CY2,CZ2,UNIT C DO 1 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)) 1 CONTINUE RETURN END C SUBROUTINE VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,N) !DEC$REAL:8 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 C CALL FGAUSS(FG,2*N,N,FFG,N) C DO 10 I=1,N VELX=DSQRT((FFG(I)*ZARG)**2+VELC) VELY=FG(I)*ZARG VELZ=FG(I+N)*ZARG C 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 10 CONTINUE RETURN END C SUBROUTINE VELOC(E,COSX,COSY,COSZ,SINE) C C FETCH A NEW VELOCITY FROM A MAXWELLIAN FLUX AT A SURFACE C !DEC$REAL:8 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 C IF (INIV1.EQ.0) CALL FGAUSS(FG,INIV1,64,FFG,INIV3) C VELX=FFG(INIV3)*ZARG VELY=FG(INIV1)*ZARG VELZ=FG(INIV1-1)*ZARG C SHEATH CONTRIBUTION IF (VELC.GT.0.) THEN VELX=DSQRT(VELC+VELX**2) ENDIF INIV1=INIV1-2 INIV3=INIV3-1 C 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 C RETURN END C SUBROUTINE FGAUSS (FG,IND,IANZ,FFG,IND2) !DEC$REAL:8 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 INTEGER IND,IND2,IANZ,JJ INTEGER*4 ISEED REAL*8 PI2,ZZ,ZSIN,ZCOS,AR,ZT REAL*8 FG(1),FFG(1) DATA PI2/6.28318530717598D0/ IND=IANZ+IANZ C CDIR$ IVDEP DO 1 JJ=1,IANZ C 1. COMPUTE THE SINE AND COSINE OF 2*PI*RAN(1) C CC ZZ=PI2*RANF() CC ZZ=PI2*DRAND48() ZZ=PI2*DBLE(RAN(ISEED)) ZSIN=DSIN(ZZ) ZCOS=DCOS(ZZ) C CC AR=DLOG(RANF()) CC AR=DLOG(DRAND48()) AR=DLOG(DBLE(RAN(ISEED))) ZT=DSQRT(-1.0D0*(AR+AR)) FG(JJ+IANZ)=ZT*ZSIN FG(JJ)=ZT*ZCOS 1 CONTINUE C C RETURN IANZ RANDOM NUMBERS FROM A GAUSSIAN FLUX IN THE ARRAY FFG C IND2=IANZ DO 2 JJ=1,IANZ CC AR=DLOG(RANF()) CC AR=DLOG(DRAND48()) AR=DLOG(DBLE(RAN(ISEED))) 2 FFG(JJ)=DSQRT(-1.D0*(AR+AR)) RETURN END C SUBROUTINE ENERGV(FE,E,COSX,COSY,COSZ,SINE,N) !DEC$REAL:8 C C FETCH A NEW ENERGY FROM A MAXWELLIAN FLUX AT A SURFACE C 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 C CALL EMAXW(FE,N) C DO 10 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 10 CONTINUE CC WRITE(6,*) (E(I),I=1,N),(COSX(I),I=1,N) RETURN END C SUBROUTINE ENERG(E,COSX,COSY,COSZ,SINE) !DEC$REAL:8 C C FETCH A NEW ENERGY FROM A MAXWELLIAN FLUX AT A SURFACE C REAL*8 FE(16) REAL*8 M1,COSX,SINE,COSY,COSZ,E,EMT REAL*8 TI,SHEATH,CALFA COMMON/B/ TI,SHEATH,CALFA C CALL EMAXW(FE,16) C 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 10 CONTINUE CC WRITE(6,*) E,COSX RETURN END C SUBROUTINE EMAXW (FE,NUMB) !DEC$REAL:8 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 INTEGER NUMB,I INTEGER*4 ISEED REAL*8 FE(1) REAL*8 PI,AR1,AR2 DATA PI/3.14159265358979D0/ C CDIR$ IVDEP DO 1 I=1,NUMB CC AR1=DLOG(RANF()) CC AR1=DLOG(DRAND48()) AR1=DLOG(DBLE(RAN(ISEED))) CC AR2=DLOG(RANF())*(DCOS(PI*0.5*RANF()))**2 CC AR2=DLOG(DRAND48())*(DCOS(PI*0.5*DRAND48()))**2 AR2=DLOG(DBLE(RAN(ISEED)))*(DCOS(PI*0.5*DBLE(RAN(ISEED))))**2 FE(I)=DSQRT(-1.D0*(AR1+AR2)) 1 CONTINUE RETURN END C FUNCTION CVMGT(A, B, C) !DEC$REAL:8 REAL*8 A,B LOGICAL C CVMGT = B IF ( C ) CVMGT = A RETURN END C SUBROUTINE SCOPY(IM,A,INCA,B,INCB) 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 10 J = 1,IM B(JB) = A(JA) JA = JA + INCA JB = JB + INCB 10 CONTINUE RETURN END C FUNCTION ILLZ(N,A,K) !DEC$REAL:8 LOGICAL A(*) INTEGER K,L,N,I INTEGER*4 ILLZ IF(K.GT.0) THEN L=N+1 DO 100 I=N,1,-1 100 IF(A(I)) L=I ELSE L=0 DO 200 I=1,N 200 IF(A(I)) L=I L=N+1-L ENDIF ILLZ=L-1 RETURN END C FUNCTION ISRCHFGE(N,ARRAY,INC,TARGT) !DEC$REAL:8 INTEGER I,N,J,INC REAL*8 ARRAY(N) REAL*8 TARGT J=1 IF(INC.LT.0) J=N*(-INC) DO 100 I=1,N IF(ARRAY(J).GE.TARGT) GO TO 200 J=J+INC 100 CONTINUE 200 ISRCHFGE=I RETURN END C FUNCTION ISRCHFGT(N,ARRAY,INC,TARGT) !DEC$REAL:8 INTEGER I,N,J,INC REAL*8 ARRAY(N),TARGT C WRITE(*,*)targt J=1 IF(INC.LT.0) J=N*(-INC) DO 100 I=1,N IF(ARRAY(J).GT.TARGT) GO TO 200 J=J+INC 100 CONTINUE 200 ISRCHFGT=I RETURN END C FUNCTION ISRCHEQ(N,ARRAY,INC,TARGT) !DEC$REAL:8 INTEGER I,N,J,INC REAL*8 ARRAY(N),TARGT C WRITE(*,*)targt J=1 IF(INC.LT.0) J=N*(-INC) DO 100 I=1,N IF(ARRAY(J).EQ.TARGT) GO TO 200 J=J+INC 100 CONTINUE 200 ISRCHEQ=I RETURN END C SUBROUTINE ENERGGAUSS(ISEED2,Esig,Epar,E0) !DEC$REAL:8 INTEGER*4 ISEED2 REAL*8 E0,Esig,Epar REAL*8 p1,p2,p3,pi DATA pi/3.14159265358979D0/ p1 = Esig*DSQRT(-2.D0*DLOG(1.D0-DBLE(RAN(ISEED2)))) p2 = 2.D0*pi*DBLE(RAN(ISEED2)) p3 = p1*DCOS(p2) Epar= E0-p3 C WRITE(*,*)E0,Esig,Epar C WRITE(31,100)E0,Esig,Epar C100 FORMAT(1x,F12.5,2x,F12.5,2x,F12.5) RETURN END C SUBROUTINE ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, + CALFA,SALFA,BW) !DEC$REAL:8 INTEGER*4 ISEED3 REAL*8 ALPHA,ALPHASIG,ALPHApar REAL*8 BW,ALFA,CALFA,SALFA REAL*8 p1,p2,p3,pi DATA pi/3.14159265358979D0/ p1 = ALPHASIG*DSQRT(-2.D0*DLOG(1.D0-DBLE(RAN(ISEED3)))) p2 = 2.D0*pi*DBLE(RAN(ISEED3)) 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) C WRITE(*,*)ALPHA,ALPHASIG,ALPHApar C WRITE(31,100)ALPHA,ALPHASIG,ALPHApar C100 FORMAT(1x,F12.5,2x,F12.5,2x,F12.5) RETURN END C 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