C 01/07/85 14:46:10 01/07/85 14:32:25 $380300.MJS SUBROUTINE DLDCUR(DN,W,U,A,Z,MODE,NCHN,LDIST,DIST) C C THIS SUBROUTINE COMPUTES THE RESPONSE OF THE D-CUP. THE INPUT C VARIABLES ARE C DN NUMBER DENSITY IN NUMBER PER CC C W THERMAL SPEED IN KM/SEC C U PLASMA BULK VELOCITY IN CUP COORDINATES C A ATOMIC MASS OF IONS IN AMU C Z CHARGE STATE OF IONS PROTON CHARGES C NCHN LIMITS ON THE CHANNEL NUMBER C THE OUTPUT VARIABLE IS C DIST THE REDUCED DISTRIBUTION FUNCTION C c REAL*8 $$$DT(6) /'$$DATE$$','01/07/85','14:32:25', c | ' DLDCUR',' FORTRAN','MJS '/ LOGICAL LDIST INTEGER*4 NSTEP(129),MODE,NCN,N INTEGER*4 NC1,NC2,NC3,NCHN(2),NCHAN,JJ,J,NM,NFIRST,NL,I C C THE INTEGER CONSTANTS HAVE THE FOLLOWING MEANING C NSTEP(I) NUMBER OF STEPS IN THE NUMERICAL INTEGRATION C BETWEEN THE THESHOLD VOLTAGES FOR CHANNEL I AND I+1 C NFIRST(I) NUMBER OF THE FIRST ELEMENT OF TVZ INCLUDED C IN THE INTEGRATION OF THE NEAR TAIL FOR CHANNEL I C REAL*4 AAAVZ,AAVZ,AVZ,AAADVZ,AADVZ,ADVZ,DVOLT,TVZ,TDVZ,DVOLTL REAL*4 AA,S,CC C REAL*4 F(3)/.6488,.1475,.07844/,D(3)/.146963,.904687,.993027/ REAL*4 F(3),D(3) REAL*8 DN,W,U(3),A,Z,DIST(128),X0,X,CONST,A0,T0 REAL*8 DVZ,G,ZZ,SS,Z0,Z1,Z10 REAL*8 ECHRGE,PI12,VT2,X2 REAL*8 GAMMA,KAPPA,CURENT,VT,G12,G32,KAPPA2,VZ,CRR(129),CRNT(129) C C THE VELOCITY SPACE INTEGRATION IS DONE IN CYLINDRICAL COORDINATES, C WITH THE CUP NORMAL POINTING IN THE Z- DIRECTION. THE C INTEGRATION OVER ANGLE IN THE VX-VY PLANE RESULTS IN THE C BESSEL FUNCTION I0; IN ORDER TO PERFORM THE NEXT INTEGRATION C ANALYTICALLY, IT IS NECESSARY TO APPROXIMATE THE BESSEL C FUNCTION BY A SUM OF EXPONENTIALS. D AND F ARE THE PARAMETERS C IN THIS APPROXIMATION. THE CARD WHICH IS COMMENTED OUT CONTAINS C ALTERNATIVE VALUES FOR THESE PARAMETERS. C C SOME OTHER REAL VARIABLES ARE: C A0 NOMINAL COLLECTING AREA OF DETECTOR C T0 NOMINAL GRID TRANSPARENCY AT NORMAL INCIDENCE C ECHRGE PROTON CHARGE IN COULOMBS C PI12 SQUARE ROOT OF PI C PRESET VALUE FOR MAXWELLIAN CUTOFF FLAG C C1,C2,C3 PARAMETERS FOR TRAPEZOIDAL APPROXIMATION TO C SENSITIVE AREA C COMMON/TRNPRD/AA(100,129,2),CC(100,129,2),S(100,129) COMMON/VZ/AAAVZ(22,5),AAVZ(4,40),AVZ(2,84),AAADVZ(22,5), , AADVZ(4,40),ADVZ(2,84),DVOLT(128),NFIRST(129),MFIRST(129) COMMON/TAIL/TVZ(241),TDVZ(241) COMMON/LCOMP/DVOLTL(16),NM(129),NL(17) data NSTEP /22,16,10,6,41*4,84*2/ data F /.6876,.2044,.09153/,D /-.22131,.851562,.990971/ data A0/82.7/,T0/.6801/ data Z0/.12D0/,Z1/1.82D0/,Z10/1.70D0/ data ECHRGE/1.6D-19/,PI12/1.7725D0/ C C DEFINE CONST AND LIMITS FOR LOOP OVER MODULATOR STEPS C CONST=DN*2.*Z*ECHRGE*A0*T0/W**3/PI12*1.D20 NC1=NCHN(1) NC2=NCHN(2) NC3=NCHN(2)+1 VT2=U(1)**2+U(2)**2 VT=DSQRT(VT2) C C ENTER OUTER LOOP OVER MODULATOR STEPS C DO 5000 NCN=NC1,NC3 NCHAN=NCN IF(MODE .EQ. 1) NCHAN=1+8*(NCN-1) C C INITIALIZE CURRENT C CURENT=0.D0 DO 2600 JJ=1,100 C C LAST INTEGRATION LOOP NOT NECCESSAY; SKIP LOOP IF APPROPRIATE C if ( ncn .ne. nc3) go to 8 if (mode -2) 4,5,8 4 if (jj .gt. nl(ncn)) go to 2600 go to 8 5 if (jj .gt. nm(ncn)) go to 2600 8 continue c IF (NCN .EQ. NC3 .AND. JJ .GT. NL(NCN) .AND. MODE .EQ. 1) c + GO TO 2600 c IF (NCN .EQ. NC3 .AND. JJ .GT. NM(NCN) .AND. MODE .EQ. 2) c + GO TO 2600 C C INTIALIZE G AND CHOOSE PROPER VALUE OF VZ,DVZ C G=0.D0 IF(NCHAN .GE. 6 .AND. NCHAN .LE. 45)GO TO 100 IF(NCHAN .GE. 46)GO TO 200 IF(JJ .GT. NSTEP(NCHAN)) GO TO 10 VZ=DBLE(AAAVZ(JJ,NCHAN))*DSQRT(Z/A) DVZ=DBLE(AAADVZ(JJ,NCHAN))*DSQRT(Z/A) GO TO 1000 10 N=JJ-NSTEP(NCHAN)+NFIRST(NCHAN)-1 VZ=DBLE(TVZ(N))*DSQRT(Z/A) DVZ=DBLE(TDVZ(N))*DSQRT(Z/A) GO TO 1000 100 IF(JJ .GT. 4) GO TO 110 VZ=DBLE(AAVZ(JJ,NCHAN-5))*DSQRT(Z/A) DVZ=DBLE(AADVZ(JJ,NCHAN-5))*DSQRT(Z/A) GO TO 1000 110 N=JJ-5+NFIRST(NCHAN) VZ=DBLE(TVZ(N))*DSQRT(Z/A) DVZ=DBLE(TDVZ(N))*DSQRT(Z/A) GO TO 1000 200 IF(JJ .GT. 2) GO TO 210 VZ=DBLE(AVZ(JJ,NCHAN-45))*DSQRT(Z/A) DVZ=DBLE(ADVZ(JJ,NCHAN-45))*DSQRT(Z/A) GO TO 1000 210 N=JJ-3+NFIRST(NCHAN) VZ=DBLE(TVZ(N))*DSQRT(Z/A) DVZ=DBLE(TDVZ(N))*DSQRT(Z/A) 1000 CONTINUE X0=((VZ-U(3))/W)**2 IF(X0 .GT. 16.0) GO TO 2650 X0=X0+VT2/W**2 C C EVALUATE INTEGRATION OVER VT ANALYTICALLY; C SUM OVER TERMS IN GRID TRANSPARENCY APPROXIMATION C DO 2500 I=1,2 GAMMA=((VZ/W)**2+DBLE(AA(JJ,NCHAN,I)))/S(JJ,NCHAN)**2 G12=DSQRT(GAMMA) C C SUM OVER TERMS IN APPROXIMATION TO BESSEL FUNCTION C DO 2500 J=1,3 KAPPA=DBLE(D(J)/S(JJ,NCHAN))*VZ*VT/W**2/G12 KAPPA2=KAPPA**2 X=X0-KAPPA2 IF(X .GT. 10)GO TO 2500 G32=G12**3 ZZ=2*G32*Z10 SS=(PI12*(KAPPA2-G12*KAPPA*Z0+.5)*DERF(G12*Z0-KAPPA)- - PI12*(KAPPA2-G12*KAPPA*Z1+.5)*DERF(G12*Z1-KAPPA))/ZZ+ + PI12*KAPPA*DERF(KAPPA)/2.D0/GAMMA X2=(G12*Z0-KAPPA)**2 IF(X2 .LT. 10)SS=SS-KAPPA*DEXP(-X2)/ZZ X2=(G12*Z1-KAPPA)**2 IF(X2 .LT. 10)SS=SS+KAPPA*DEXP(-X2)/ZZ IF(KAPPA2 .LT. 10)SS=SS+DEXP(-KAPPA2)/2.D0/GAMMA C C ADD TO SUM FOR INTEGRAND (MAXWELLIAN FACTOR NOT INCLUDED) C G=G+SS*DBLE(F(J)*CC(JJ,NCHAN,I))*DEXP(-X) 2500 CONTINUE C C ADD TO SUM FOR NUMERICAL INTEGRATION OVER VZ C CURENT=CURENT+G*(VZ/S(JJ,NCHAN))**2*DVZ*VZ C C FILL ARRAY FOR UPPER MODULATOR STEP INTEGRATION, IF APPROPRIATE C 2650 CONTINUE if (mode - 2) 2591,2592,2595 2591 if (jj .eq. nl(ncn)) crr(ncn) = curent*const go to 2595 2592 if ( jj .eq. nm(ncn)) crr(ncn) = curent*const 2595 continue c IF (JJ .EQ. NL(NCN) .AND. MODE .EQ. 1) CRR(NCN)=CURENT*CONST c IF (JJ .EQ. NM(NCN) .AND. MODE .EQ. 2) CRR(NCN)=CURENT*CONST 2600 CONTINUE C C FILL ARRAY FOR LOWER MODULATOR STEP INTEGRATION C CRNT(NCN)=CURENT*CONST 5000 CONTINUE C C COMPUTE CURRENTS BY TAKING DIFFERENCE BETWEEN ADJACENT MODULATOR C STEPS; IF LDIST= .TRUE., DIVIDE BY CHANNEL WIDTH C DO 5110 JJ=NC1,NC2 IF(MODE .EQ. 2) GO TO 5100 C C L-MODE; LDIST = .TRUE. C IF(LDIST)DIST(JJ)=(CRNT(JJ)-CRR(JJ+1))/DBLE(DVOLTL(JJ)) GO TO 5110 C C L-MODE; LDIST = .TRUE. C 5100 IF(LDIST)DIST(JJ)=(CRNT(JJ)-CRR(JJ+1))/DBLE(DVOLT(JJ)) C C L- OF M-MODE; LDIST = .FALSE. C 5110 IF(.NOT. LDIST)DIST(JJ)=(CRNT(JJ)-CRR(JJ+1)) RETURN END