diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 1bc3684..386742c 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -2582,67 +2582,11 @@ C C C BACKSCATTERING C - IF(IB.EQ.0) GO TO 1512 - WRITE(21,1527) - 1527 FORMAT(1H1,//5X,'BACKSCATTERING OF PROJECTILES'/) - BI=DBLE(IB) - BIL=DBLE(IBL) - RN=BI/HN - RE=EB/(HN*EMV) - EMEANR=RE/RN - EMEAN=EB/BI - AVEB=EMEAN - IF (equal(BI,1.0d0))GO TO 1506 - AVNLB=ENUCLB/BI - VANLB=ENL2B/BI-AVNLB*AVNLB - SIGNLB=DSQRT(VANLB) - DFINLB=SIGNLB/BI - AVILB=EINELB/BI - VAILB=EIL2B/BI-AVILB*AVILB - SIGILB=DSQRT(VAILB) - DFIILB=SIGILB/BI - 1506 WRITE(21,1508) RN,RE,EMEANR,EMEAN - 1508 FORMAT(/5X,'PART.REFL.COEF.=',1PE11.4,' ENERGY REFL.COEF.=' - & ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' - & ,1E11.4) - IF(IB.EQ.0) GO TO 1512 - CALL MOMENT(EB1B,EB2B,EB3B,EB4B,EB5B,EB6B,EB,EB2SUM,EB3SUM,EB4SUM - & ,EB5SUM,EB6SUM,BI) - CALL MOMENT(EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL,EB1SUL,EB2SUL - & ,EB3SUL,EB4SUL,EB5SUL,EB6SUL,BIL) - CALL MOMENTS(FIB0,SEB,THB,FOB,FIB,SIB,SIGMAB,DFIB0,DSEB,DTHB,EB - & ,EB2SUM,EB3SUM,EB4SUM,EB5SUM,EB6SUM,BI) - CALL MOMENT(PL1S,PL2S,PL3S,PL4S,PL5S,PL6S,PLSB,PL2SB,PL3SB,PL4SB - & ,PL5SB,PL6SB,BI) - CALL MOMENTS(FIPB0,SEPB,THPB,FOPB,FIPB,SIPB,SIGMPB,DFIPB0,DSEPB - & ,DTHPB, PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB,BI) - WRITE(21,7117) - WRITE(21,7241) FIB0,SEB,THB,FOB,SIGMAB,DFIB0,DSEB,DTHB - 7241 FORMAT(1X,' ENERGY',5X,1P1E12.4,7E14.4) - WRITE(21,7242) FIPB0,SEPB,THPB,FOPB,SIGMPB,DFIPB0,DSEPB,DTHPB - 7242 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) - WRITE(21,7237) AVNLB,VANLB,SIGNLB,DFINLB - WRITE(21,7238) AVILB,VAILB,SIGILB,DFIILB - WRITE(21,7118) - WRITE(21,1541) EB1B,EB2B,EB3B,EB4B,EB5B,EB6B - 1541 FORMAT(1X,' ENERGY',5X,1P1E12.4,5E14.4) - WRITE(21,1543) EB1BL,EB2BL,EB3BL,EB4BL,EB5BL,EB6BL - 1543 FORMAT(1X,' LOGENERGY',5X,1P1E12.4,5E14.4) - WRITE(21,1545) PL1S,PL2S,PL3S,PL4S,PL5S,PL6S - 1545 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) - DO I=1,20 - AI(I)=0.05D0*DBLE(I) - ENDDO - IF(IB.EQ.0) GO TO 1512 - WRITE(21,1514) - 1514 FORMAT(//5X,'POLAR ANGULAR DISTRIBUTION OF BACKSCATTERED ', - & 'PROJECTILES'//) - DO I=1,20 - RKADB(I)=DBLE(KADB(I))*20.D0/DBLE(IB) - ENDDO - WRITE(21,1518)(AI(I),I=1,20),(KADB(I),I=1,20),(RKADB(I),I=1,20) - 1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) - 1512 CONTINUE + call writeBackscatterSummary(IB,IBL,HN,EB,EMV,ENUCLB, + & ENL2B,EINELB,EIL2B,EB2SUM,EB3SUM,EB4SUM,EB5SUM, + & EB6SUM,EB1SUL,EB2SUL,EB3SUL,EB4SUL,EB5SUL, + & EB6SUL,PLSB,PL2SB,PL3SB,PL4SB,PL5SB,PL6SB, + & KADB,AI,RKADB,FIB0,SIGMAB) C C TRANSMISSION C @@ -2667,6 +2611,9 @@ C 1522 FORMAT(//5X,'PART.TRANSM.COEF.=',1PE11.4,' ENERGY TRANSM.COEF.=' & ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' & ,1E11.4) + 7241 FORMAT(1X,' ENERGY',5X,1P1E12.4,7E14.4) + 7242 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + 1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) CALL MOMENTS(FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT,ET & ,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM,TIT) CALL MOMENTS(FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT,DFIPT0,DSEPT @@ -3650,6 +3597,149 @@ C======================================================================= C OUTPUT WRITER SUBROUTINES C======================================================================= +C======================================================================= +C WRITE BACKSCATTERING SUMMARY +C======================================================================= + subroutine writeBackscatterSummary(nBackscatter,nBackLayer, + & nPrimary,backEnergySum,energyUnit,elasticLossSum, + & elasticLossSquare,inElasticLossSum,inElasticLossSquare, + & energyMoment2,energyMoment3,energyMoment4,energyMoment5, + & energyMoment6,logEnergyMoment1,logEnergyMoment2, + & logEnergyMoment3,logEnergyMoment4,logEnergyMoment5, + & logEnergyMoment6,pathMoment1,pathMoment2,pathMoment3, + & pathMoment4,pathMoment5,pathMoment6,angleBins,angleX, + & angleDistribution,meanBackEnergy,backEnergySigma) + implicit none + integer nBackscatter,nBackLayer + integer angleBins(20) + integer i + real*8 nPrimary,backEnergySum,energyUnit + real*8 elasticLossSum,elasticLossSquare + real*8 inElasticLossSum,inElasticLossSquare + real*8 energyMoment2,energyMoment3,energyMoment4 + real*8 energyMoment5,energyMoment6 + real*8 logEnergyMoment1,logEnergyMoment2,logEnergyMoment3 + real*8 logEnergyMoment4,logEnergyMoment5,logEnergyMoment6 + real*8 pathMoment1,pathMoment2,pathMoment3 + real*8 pathMoment4,pathMoment5,pathMoment6 + real*8 angleX(20),angleDistribution(20) + real*8 meanBackEnergy,backEnergySigma + real*8 bi,bil,particleReflection,energyReflection + real*8 relativeMeanEnergy,meanEnergy + real*8 meanElasticLoss,varElasticLoss,sigmaElasticLoss + real*8 errorElasticLoss + real*8 meanInElasticLoss,varInElasticLoss,sigmaInElasticLoss + real*8 errorInElasticLoss + real*8 energyVariance,energySkew,energyKurtosis + real*8 energyFifth,energySixth + real*8 energyError1,energyError2,energyError3 + real*8 pathMean,pathVariance,pathSkew,pathKurtosis + real*8 pathFifth,pathSixth,pathSigma + real*8 pathError1,pathError2,pathError3 + real*8 energyMoment1Out,energyMoment2Out,energyMoment3Out + real*8 energyMoment4Out,energyMoment5Out,energyMoment6Out + real*8 logMoment1Out,logMoment2Out,logMoment3Out + real*8 logMoment4Out,logMoment5Out,logMoment6Out + real*8 pathMoment1Out,pathMoment2Out,pathMoment3Out + real*8 pathMoment4Out,pathMoment5Out,pathMoment6Out + logical equal + + if(nBackscatter.eq.0) return + write(21,1527) + 1527 format(1H1,//5X,'BACKSCATTERING OF PROJECTILES'/) + bi=dble(nBackscatter) + bil=dble(nBackLayer) + particleReflection=bi/nPrimary + energyReflection=backEnergySum/(nPrimary*energyUnit) + relativeMeanEnergy=energyReflection/particleReflection + meanEnergy=backEnergySum/bi + meanBackEnergy=meanEnergy + if (equal(bi,1.0d0)) go to 1506 + meanElasticLoss=elasticLossSum/bi + varElasticLoss=elasticLossSquare/bi + & -meanElasticLoss*meanElasticLoss + sigmaElasticLoss=dsqrt(varElasticLoss) + errorElasticLoss=sigmaElasticLoss/bi + meanInElasticLoss=inElasticLossSum/bi + varInElasticLoss=inElasticLossSquare/bi + & -meanInElasticLoss*meanInElasticLoss + sigmaInElasticLoss=dsqrt(varInElasticLoss) + errorInElasticLoss=sigmaInElasticLoss/bi + 1506 write(21,1508) particleReflection,energyReflection, + & relativeMeanEnergy,meanEnergy + 1508 format(/5X,'PART.REFL.COEF.=',1PE11.4, + & ' ENERGY REFL.COEF.=',1E11.4, + & ' REL.MEAN ENERGY =',1E11.4, + & ' MEAN ENERGY =',1E11.4) + call MOMENT(energyMoment1Out,energyMoment2Out,energyMoment3Out, + & energyMoment4Out,energyMoment5Out,energyMoment6Out, + & backEnergySum,energyMoment2,energyMoment3,energyMoment4, + & energyMoment5,energyMoment6,bi) + call MOMENT(logMoment1Out,logMoment2Out,logMoment3Out, + & logMoment4Out,logMoment5Out,logMoment6Out, + & logEnergyMoment1,logEnergyMoment2,logEnergyMoment3, + & logEnergyMoment4,logEnergyMoment5,logEnergyMoment6,bil) + call MOMENTS(meanBackEnergy,energyVariance,energySkew, + & energyKurtosis,energyFifth,energySixth,backEnergySigma, + & energyError1,energyError2,energyError3,backEnergySum, + & energyMoment2,energyMoment3,energyMoment4,energyMoment5, + & energyMoment6,bi) + call MOMENT(pathMoment1Out,pathMoment2Out,pathMoment3Out, + & pathMoment4Out,pathMoment5Out,pathMoment6Out,pathMoment1, + & pathMoment2,pathMoment3,pathMoment4,pathMoment5, + & pathMoment6,bi) + call MOMENTS(pathMean,pathVariance,pathSkew,pathKurtosis, + & pathFifth,pathSixth,pathSigma,pathError1,pathError2, + & pathError3,pathMoment1,pathMoment2,pathMoment3,pathMoment4, + & pathMoment5,pathMoment6,bi) + write(21,7117) + write(21,7241) meanBackEnergy,energyVariance,energySkew, + & energyKurtosis,backEnergySigma,energyError1,energyError2, + & energyError3 + 7241 format(1X,' ENERGY',5X,1P1E12.4,7E14.4) + write(21,7242) pathMean,pathVariance,pathSkew,pathKurtosis, + & pathSigma,pathError1,pathError2,pathError3 + 7242 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + write(21,7237) meanElasticLoss,varElasticLoss,sigmaElasticLoss, + & errorElasticLoss + write(21,7238) meanInElasticLoss,varInElasticLoss, + & sigmaInElasticLoss,errorInElasticLoss + write(21,7118) + write(21,1541) energyMoment1Out,energyMoment2Out, + & energyMoment3Out,energyMoment4Out,energyMoment5Out, + & energyMoment6Out + 1541 format(1X,' ENERGY',5X,1P1E12.4,5E14.4) + write(21,1543) logMoment1Out,logMoment2Out,logMoment3Out, + & logMoment4Out,logMoment5Out,logMoment6Out + 1543 format(1X,' LOGENERGY',5X,1P1E12.4,5E14.4) + write(21,1545) pathMoment1Out,pathMoment2Out,pathMoment3Out, + & pathMoment4Out,pathMoment5Out,pathMoment6Out + 1545 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + do i=1,20 + angleX(i)=0.05d0*dble(i) + enddo + write(21,1514) + 1514 format(//5X,'POLAR ANGULAR DISTRIBUTION OF BACKSCATTERED ', + & 'PROJECTILES'//) + do i=1,20 + angleDistribution(i)=dble(angleBins(i))*20.d0 + & /dble(nBackscatter) + enddo + write(21,1518)(angleX(i),i=1,20),(angleBins(i),i=1,20), + & (angleDistribution(i),i=1,20) + 1518 format(5X,20F6.2//,5X,20I6/5X,20F6.3) + return + 7117 format(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ', + & 4X,' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ', + & 3X,' ERROR 2.M ',3X,' ERROR 3.M ') + 7118 format(/20X,' 1.MOMENT ',4X,' 2.MOMENT ',4X,' 3.MOMENT ', + & 4X,' 4.MOMENT ',4X,' 5.MOMENT ',4X,' 6.MOMENT ') + 7237 format(1X,'ELASTIC LOSS',5X,1P1E12.4,1E14.4,28X, + & 2E14.4) + 7238 format(1X,' INEL. LOSS',5X,1P1E12.4,1E14.4,28X, + & 2E14.4) + end + C======================================================================= C WRITE IMPLANTATION PROFILE AND RECOIL DEPTH TABLES C=======================================================================