diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 386742c..2f9cbd0 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -2590,48 +2590,11 @@ C C C TRANSMISSION C - IF(IT.EQ.0) GO TO 1524 - WRITE(21,1529) - 1529 FORMAT(///5X,' TRANSMISSION OF PROJECTILES'/) - TIT=DBLE(IT) - TN=TIT/HN - TE=ET/(HN*E0) - TMEANR=TE/TN - EMEANT=TMEANR*E0 - IF (equal(TIT,1.0D0)) GO TO 1520 - AVNLT=ENUCLT/TIT - VANLT=ENL2T/TIT-AVNLT*AVNLT - SIGNLT=DSQRT(VANLT) - DFINLT=SIGNLT/TIT - AVILT=EINELT/TIT - VAILT=EIL2T/TIT-AVILT*AVILT - SIGILT=DSQRT(VAILT) - DFIILT=SIGILT/TIT - 1520 WRITE(21,1522) TN,TE,TMEANR,EMEANT - 1522 FORMAT(//5X,'PART.TRANSM.COEF.=',1PE11.4,' ENERGY TRANSM.COEF.=' - & ,1E11.4,' REL.MEAN ENERGY =',1E11.4,' MEAN ENERGY =' - & ,1E11.4) - 7241 FORMAT(1X,' ENERGY',5X,1P1E12.4,7E14.4) - 7242 FORMAT(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + call writeTransmissionSummary(IT,HN,ET,E0,ENUCLT,ENL2T, + & EINELT,EIL2T,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM, + & PLST,PL2ST,PL3ST,PL4ST,PL5ST,PL6ST,KADT,AI,RKADT, + & FIT0,SIGMAT) 1518 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) - CALL MOMENTS(FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT,ET - & ,ET2SUM,ET3SUM,ET4SUM,ET5SUM,ET6SUM,TIT) - CALL MOMENTS(FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT,DFIPT0,DSEPT - & ,DTHPT, PLST,PL2ST,PL3ST,PL4ST,PL5ST,PL6ST,TIT) - WRITE(21,7117) - WRITE(21,7241) FIT0,SET,THT,FOT,SIGMAT,DFIT0,DSET,DTHT - WRITE(21,7242) FIPT0,SEPT,THPT,FOPT,SIGMPT,DFIPT0,DSEPT,DTHPT - WRITE(21,7237) AVNLT,VANLT,SIGNLT,DFINLT - WRITE(21,7238) AVILT,VAILT,SIGILT,DFIILT - WRITE(21,1526) - 1526 FORMAT(//5X,'POLAR ANGULAR DISTRIBUTION OF TRANSMITTED ', - & 'PARTICLES'//) - DO I=1,20 - RKADT(I)=DBLE(KADT(I))*20.D0/DBLE(IT) - ENDDO - WRITE(21,1530) (AI(I),I=1,20),(KADT(I),I=1,20),(RKADT(I),I=1,20) - 1530 FORMAT(5X,20F6.2//,5X,20I6/5X,20F6.3) - 1524 CONTINUE C C BACKWARD SPUTTERING : YIELDS AND ENERGIES C @@ -3597,6 +3560,117 @@ C======================================================================= C OUTPUT WRITER SUBROUTINES C======================================================================= +C======================================================================= +C WRITE TRANSMISSION SUMMARY +C======================================================================= + subroutine writeTransmissionSummary(nTrans,nPrimary,transEnergy, + & initialEnergy,elasticLossSum,elasticLossSquare, + & inElasticLossSum,inElasticLossSquare,energyMoment2, + & energyMoment3,energyMoment4,energyMoment5,energyMoment6, + & pathMoment1,pathMoment2,pathMoment3,pathMoment4, + & pathMoment5,pathMoment6,angleBins,angleX, + & angleDistribution,meanTransEnergy,transEnergySigma) + implicit none + integer nTrans + integer angleBins(20) + integer i + real*8 nPrimary,transEnergy,initialEnergy + real*8 elasticLossSum,elasticLossSquare + real*8 inElasticLossSum,inElasticLossSquare + real*8 energyMoment2,energyMoment3,energyMoment4 + real*8 energyMoment5,energyMoment6 + real*8 pathMoment1,pathMoment2,pathMoment3 + real*8 pathMoment4,pathMoment5,pathMoment6 + real*8 angleX(20),angleDistribution(20) + real*8 meanTransEnergy,transEnergySigma + real*8 transCount,particleTransmission,energyTransmission + real*8 relativeMeanEnergy,meanEnergy + real*8 meanElasticLoss,varElasticLoss,sigmaElasticLoss + real*8 errorElasticLoss + real*8 meanInElasticLoss,varInElasticLoss,sigmaInElasticLoss + real*8 errorInElasticLoss + real*8 energyVariance,energySkew,energyKurtosis + real*8 energyFifth,energySixth + real*8 energyError1,energyError2,energyError3 + real*8 pathMean,pathVariance,pathSkew,pathKurtosis + real*8 pathFifth,pathSixth,pathSigma + real*8 pathError1,pathError2,pathError3 + logical equal + + if(nTrans.eq.0) return + write(21,1529) + 1529 format(///5X,' TRANSMISSION OF PROJECTILES'/) + transCount=dble(nTrans) + particleTransmission=transCount/nPrimary + energyTransmission=transEnergy/(nPrimary*initialEnergy) + relativeMeanEnergy=energyTransmission/particleTransmission + meanEnergy=relativeMeanEnergy*initialEnergy + meanTransEnergy=meanEnergy + meanElasticLoss=0.0d0 + varElasticLoss=0.0d0 + sigmaElasticLoss=0.0d0 + errorElasticLoss=0.0d0 + meanInElasticLoss=0.0d0 + varInElasticLoss=0.0d0 + sigmaInElasticLoss=0.0d0 + errorInElasticLoss=0.0d0 + if (equal(transCount,1.0d0)) go to 1520 + meanElasticLoss=elasticLossSum/transCount + varElasticLoss=elasticLossSquare/transCount + & -meanElasticLoss*meanElasticLoss + sigmaElasticLoss=dsqrt(varElasticLoss) + errorElasticLoss=sigmaElasticLoss/transCount + meanInElasticLoss=inElasticLossSum/transCount + varInElasticLoss=inElasticLossSquare/transCount + & -meanInElasticLoss*meanInElasticLoss + sigmaInElasticLoss=dsqrt(varInElasticLoss) + errorInElasticLoss=sigmaInElasticLoss/transCount + 1520 write(21,1522) particleTransmission,energyTransmission, + & relativeMeanEnergy,meanEnergy + 1522 format(//5X,'PART.TRANSM.COEF.=',1PE11.4, + & ' ENERGY TRANSM.COEF.=',1E11.4, + & ' REL.MEAN ENERGY =',1E11.4, + & ' MEAN ENERGY =',1E11.4) + call MOMENTS(meanTransEnergy,energyVariance,energySkew, + & energyKurtosis,energyFifth,energySixth,transEnergySigma, + & energyError1,energyError2,energyError3,transEnergy, + & energyMoment2,energyMoment3,energyMoment4,energyMoment5, + & energyMoment6,transCount) + call MOMENTS(pathMean,pathVariance,pathSkew,pathKurtosis, + & pathFifth,pathSixth,pathSigma,pathError1,pathError2, + & pathError3,pathMoment1,pathMoment2,pathMoment3,pathMoment4, + & pathMoment5,pathMoment6,transCount) + write(21,7117) + write(21,7241) meanTransEnergy,energyVariance,energySkew, + & energyKurtosis,transEnergySigma,energyError1, + & energyError2,energyError3 + 7241 format(1X,' ENERGY',5X,1P1E12.4,7E14.4) + write(21,7242) pathMean,pathVariance,pathSkew,pathKurtosis, + & pathSigma,pathError1,pathError2,pathError3 + 7242 format(1X,' PATHLENGTH',5X,1P1E12.4,7E14.4) + write(21,7237) meanElasticLoss,varElasticLoss,sigmaElasticLoss, + & errorElasticLoss + write(21,7238) meanInElasticLoss,varInElasticLoss, + & sigmaInElasticLoss,errorInElasticLoss + write(21,1526) + 1526 format(//5X,'POLAR ANGULAR DISTRIBUTION OF TRANSMITTED ', + & 'PARTICLES'//) + do i=1,20 + angleDistribution(i)=dble(angleBins(i))*20.d0/dble(nTrans) + enddo + write(21,1530) (angleX(i),i=1,20),(angleBins(i),i=1,20), + & (angleDistribution(i),i=1,20) + 1530 format(5X,20F6.2//,5X,20I6/5X,20F6.3) + return + 7117 format(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ', + & 4X,' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ', + & 3X,' ERROR 2.M ',3X,' ERROR 3.M ') + 7237 format(1X,'ELASTIC LOSS',5X,1P1E12.4,1E14.4,28X, + & 2E14.4) + 7238 format(1X,' INEL. LOSS',5X,1P1E12.4,1E14.4,28X, + & 2E14.4) + end + C======================================================================= C WRITE BACKSCATTERING SUMMARY C=======================================================================