diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index f0b145f..299b598 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -206,7 +206,7 @@ C CH1,2,3,4,5 are values of A-1,2,3,4,5 of the Ziegler tables REAL*8 E2,AB,FP,AN REAL*8 Esig,Epar REAL*8 E0keV,EsigkeV -c fuer part. reflec. coeff. Berechnung +c for part. reflec. coeff. calculation REAL*8 epsilon, prcoeff,PRC(6) REAL*8 cossin,rphi,vanlb,vailb,vanlt,vailt,phip,ta,ta2 REAL*8 exir,phipr,u,rq,es,vanli,vaili @@ -490,11 +490,12 @@ C SET INTERVAL CONSTANTS FOR OUTPUT 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 + IF(E0.LT.0.) then + E0DE = 10.0D0/(DABS(E0)*DE) + else + E0DE = 100.0D0/(E0*DE) + endif + BW = 180.D0/PI DAW = BW/DA DGW = BW/DG DGIK = BW/DGI @@ -511,25 +512,26 @@ C is wrong. Calculation form layer thickness should be correct. C The rge files (stopping profiles) and summary files should be fine. C CW is depth increment and DX(K) is the thickness of layer K depth_interval_flag = 1 - DO K=1,L-1 - write(*,*) DX(K),CW - IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN +C DO K=1,L-1 +C IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN C If the thickness of layer K is not an integer multiples of depth increment - depth_interval_flag = 0 - write (*,*) 'I am here' - GO TO 44 - ENDIF - ENDDO - 44 CONTINUE +C depth_interval_flag = 0 +C GO TO 44 +C ENDIF +C ENDDO +C 44 CONTINUE JT(1) = 0 JT(2) = NJ(1) LJ= NJ(1)+NJ(2) XX(1)=DX(1) -C Loop over all defined layers NJJ=0 +C Loop over all defined layers DO I=1,L + if(.not.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) then + depth_interval_flag = 0 + endif if (I.ge.2) XX(I)=XX(I-1)+DX(I) if (I.ge.3) then JT(I)=JT(I-1)+NJ(I-1) @@ -668,6 +670,7 @@ C SET CONSTANT DISTANCES XC = CVMGT( X0, -SU, X0.GT.0.D0) RT = TT-RD + IF(E0.GE.0.D0) GO TO 51 C C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION @@ -693,21 +696,16 @@ 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 IV=1,NUM EMX = EMX+E(IV) - ENDDO - DO iv=1,num ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) +C me(ne) is undefined here! + write (*,*) 'me=',me(ne) me(ne) = me(ne)+1 ENDDO GO TO 56 -C -C MAXWELLIAN ENERGY DISTRIBUTION -C ELSE CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM) DO IV=1,NUM @@ -823,8 +821,6 @@ C Y(IV) = LM(JL) *RA *COSY(IV) Z(IV) = LM(JL) *RA *COSZ(IV) PL(IV) = CVMGT(0.D0,LM(JL)*RA,X0.LE.0.0) - ENDDO - DO IV=1,NUM LLL(IV) = JL ENDDO C @@ -852,10 +848,7 @@ C DO IV=1,IH1 call ranlux(random, 1) JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1, - & DBLE(random(1))) - & +JT(LLL(IV)) - ENDDO - DO IV=1,IH1 + & DBLE(random(1)))+JT(LLL(IV)) EPS(IV)=E(IV)*F1(JJJ(IV)) ENDDO DO IV=1,IH1 @@ -1149,8 +1142,6 @@ C ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV)))) ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1))) ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0)) - ENDDO - DO IV=1,IH1 DEN2=DEN(IV)*DEN(IV) DEN3=DEN2*DEN(IV) EEL=EEL+DEN(IV) @@ -1281,9 +1272,6 @@ C & *CSYR(IREC1,2) ZR(IREC1,2)=ZR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2)) & *CSZR(IREC1,2) - ENDDO - - DO IREC1=1,NREC2 CXR(IREC1)=CSXR(IREC1,2) CYR(IREC1)=CSYR(IREC1,2) CZR(IREC1)=CSZR(IREC1,2) @@ -1564,8 +1552,6 @@ C ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= ICDIRJ(JJR(IREC1,2) & ,JJR(IREC1,1)) +IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT & .DI(JJR(IREC1,1)))) - ENDDO - DO IREC1=1,NREC2 TEL=TEL+TS(IREC1) TR=TR+T(IREC1) EI=EI+DEER(IREC1) @@ -1853,7 +1839,7 @@ C I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), MAXD1), 0) IRP(I1)=IRP(I1)+1 C -C Calculated the stopped particles in each layer +C Calculate the stopped particles in each layer C LowTiefe = 0.D0 UpTiefe = DX(1) @@ -2250,8 +2236,6 @@ C IF(IVMAX.LT.IVMIN) GO TO 132 DO IV=IVMAX+1,IH1 LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) - ENDDO - DO IV=IVMAX+1,IH1 X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) @@ -2272,14 +2256,14 @@ C ima = 1 1010 ima = MIN0(ima+2,5000) C I am not sure what is this file for?! - open(20,file='edist') - do ne=1,ima - write(20,1020) ne,me(ne) - enddo - 1020 format(1x,2i6) - close(20) +C open(20,file='edist') +C do ne=1,ima +C write(20,1020) ne,me(ne) +C enddo +C 1020 format(1x,2i6) +C close(20) C -C Berechnung der part. reflec. coeff. nach Thomas et al. +C Calculate the part. reflec. coeff. according to Thomas et al. C E0keV=E0/1.D3 EsigkeV=Esig/1.D3 @@ -2494,8 +2478,6 @@ C DO J=1,LJ ISPA = ISPA+IBSP(J) ESPA = ESPA+EBSP(J) - ENDDO - DO J=1,LJ ISPAT = ISPAT+ITSP(J) ESPAT = ESPAT+ETSP(J) ENDDO @@ -2515,7 +2497,9 @@ C nli(i)=IDINT(xx(i)/cw+0.01) enddo do i=1,l +C we have a problem here, what is nli(i-1) for i=1 (assumed 0??) do j=nli(i-1)+1,nli(i) + write(*,*) 'i=',i,'j=',j,nli(i-1)+1,nli(i) irpl(i)=irpl(i)+irp(j) enddo enddo @@ -3388,9 +3372,11 @@ C & ,(RKDSTJ(I,J),I=1,20) ENDDO 1800 CONTINUE -c -c The file for33 is created here -c + +C +C The file fort.33 is created here +C This is the summary file writing procedure +C 7802 FORMAT(6x,'Energy',4x,'SigmaE',5x,'Alpha',2x,'SigAlpha',4x,'ntot', & 5x,'imp',2x,'backsc',3x,'trans',3x,'tried',4x,'negE',3x, & 'range',6x,'straggeling',2x, 'Eback',7x,'sigEback',4x @@ -3403,8 +3389,6 @@ C & I0,3x)) inquire(FILE='fort.33',EXIST=FORT33) if (.not.FORT33) then open(33) -C This is where we can insert the chemical formula -C WRITE(33,7802) ('impL',k,k=1,NLayers) WRITE(33,7802) (chem(k),k=1,NLayers) else open(33,access='append') @@ -3413,8 +3397,9 @@ C WRITE(33,7802) ('impL',k,k=1,NLayers) & ,negE,FIX0,SIGMAX,FIB0,SIGMAB,FIT0,SIGMAT,epsilon,prcoeff & ,(number_in_layer(k),k=1,NLayers) CLOSE(33) -c -c End of file for33 +C +C End of file fort.33 + C C TOP AND FRONT LINES FOR MATRICES C diff --git a/trimspNL b/trimspNL index eb58964..8e64e10 100755 Binary files a/trimspNL and b/trimspNL differ