REAL*4 FUNCTION shieldose2idl(argc, argv) ! Called by IDL INCLUDE 'wrappers.inc' c INTEGER*4 argc, argv(*) ! Argc and Argv are integers j = loc(argc) ! Obtains the number of arguments (argc) ! Because argc is passed by VALUE. call shieldose2(%VAL(argv(1)),%VAL(argv(2)),%VAL(argv(3)), &%VAL(argv(4)),%VAL(argv(5)),%VAL(argv(6)),%VAL(argv(7)), &%VAL(argv(8)),%VAL(argv(9)),%VAL(argv(10)),%VAL(argv(11)), &%VAL(argv(12)),%VAL(argv(13)),%VAL(argv(14)),%VAL(argv(15)), &%VAL(argv(16)),%VAL(argv(17)),%VAL(argv(18)),%VAL(argv(19)), &%VAL(argv(20)),%VAL(argv(21)),%VAL(argv(22)),%VAL(argv(23)), &%VAL(argv(24)),%VAL(argv(25)),%VAL(argv(26)),%VAL(argv(27)), &%VAL(argv(28)),%VAL(argv(29))) shieldose2idl = 9.9 RETURN END ! PROGRAM SD2 C C SHIELDOSE-2, VERSION 2.10, 28 APR 94. C C S.M. SELTZER C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C (301) 975-5552 C C IDET = 1, AL DETECTOR C 2, GRAPHITE DETECTOR C 3, SI DETECTOR C 4, AIR DETECTOR C 5, BONE DETECTOR C 6, CALCIUM FLUORIDE DETECTOR C 7, GALLIUM ARSENIDE DETECTOR C 8, LITHIUM FLUORIDE DETECTOR C 9, SILICON DIOXIDE DETECTOR C 10, TISSUE DETECTOR C 11, WATER DETECTOR C C INUC = 1, NO NUCLEAR ATTENUATION FOR PROTONS IN AL C 2, NUCLEAR ATTENUATION, LOCAL CHARGED-SECONDARY ENERGY C DEPOSITION C 3, NUCLEAR ATTENUATION, LOCAL CHARGED-SECONDARY ENERGY C DEPOSITION, AND APPROX EXPONENTIAL DISTRIBUTION OF C NEUTRON DOSE C C IMAX = Number of shielding depth (max allowed=71) C C IUNT = the shielding depth unit C 1, Mils C 2, g/cm2 C 3, mm C C Zin = array(IMAX) of thickness in unit of IUNT C C EMINS,EMAXS= the min and max energy of solar protons spectrum [MeV] C EMINP,EMAXP= the min and max energy of trapped protons spectrum [MeV] C NPTSP= number of spectrum points which divides proton spectra for integration. C EMINE,EMAXE= the min and max energy of trapped electrons spectrum [MeV] C NPTSE=number of spectrum points which divides electron spectrum for integration. C C JSMAX= number of points in falling spectrum of solar protons C JPMAX= number of points in falling spectrum of trapped protons C JEMAX= number of points in falling spectrum of trapped electrons C C EUNIT= CONVERSION FACTOR FROM /ENERGY TO /MEV, C E.G., EUNIT = 1000 IF FLUX IS /KEV. C C DURATN= MISSION DURATION IN MULTIPLES OF UNIT TIME (s). C C INCIDENT OMNIDIRECTIONAL FLUX IN /ENERGY/CM2/UNIT TIME C (SOLAR-FLARE FLUX IN /ENERGY/CM2). C ESin: Energy array(JSMAX) of solar proton spectrum C SFLUXin: Flux array(JSMAX) for solar protons( SOLAR-FLARE FLUX IN /ENERGY/CM2) C C EPin: Energy array(JPMAX) of trapped proton spectrum C PFLUXin: Flux array(JPMAX) for trapped protons (INCIDENT OMNIDIRECTIONAL FLUX IN /ENERGY/CM2/UNIT TIME) C C EEin: Energy array(JEMAX) of trapped electron spectrum C EFLUXin: Flux array(JEMAX) for trapped electrons (INCIDENT OMNIDIRECTIONAL FLUX IN /ENERGY/CM2/UNIT TIME) C C OUTPUTS: C SolDose: Dose profile array (Imax,3) for solar protons C ProtDose: Dose profile array (Imax,3) for trapped protons C ElecDose: Dose profile array (Imax,3) for trapped electrons C BremDose: Dose profile array (Imax,3) for Brem C TotDose: Total dose profile array (Imax,3) C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) subroutine shieldose2(IDET,INUC,IMAX,IUNT,Zin, * EMINS,EMAXS,EMINP,EMAXP,NPTSP,EMINE,EMAXE,NPTSE, * JSMAX,JPMAX,JEMAX,EUNIT,DURATN, * ESin,SFLUXin,EPin,PFLUXin,EEin,EFLUXin, * SolDose,ProtDose,ElecDose,BremDose,TotDose) implicit none integer*4 MMAXPI,KMAXPI,NMAXPI,LMAXPI,IMIXI,MMAXEI,NMAXEI,LMAXSI integer*4 LMAXEI,LMAXTI,LMAXBI,IMAXI,NPTSI,JMAXI PARAMETER (MMAXPI=133,KMAXPI=30,NMAXPI=49,LMAXPI=51,IMIXI=11, 1 MMAXEI=81,NMAXEI=14,LMAXSI=33+1,LMAXEI=51,LMAXTI=37,LMAXBI=47, 2 IMAXI=71,NPTSI=1001,JMAXI=301) integer*4 NPTSPI,NPTSEI PARAMETER (NPTSPI=NPTSI,NPTSEI=NPTSI) real*8 ZCON,ZMCON PARAMETER (ZCON=0.001d0*2.540005d0*2.70d0,ZMCON=10.0d0/2.70d0) CHARACTER TAG*81,DET(IMIXI)*8, 1 VERSION*4 real*8 EP(MMAXPI),RP(MMAXPI),RPB(MMAXPI),RPC(MMAXPI), 1 RPD(MMAXPI),TEPN(KMAXPI),FEPN(KMAXPI),FEPNB(KMAXPI), 2 FEPNC(KMAXPI),FEPND(KMAXPI),TP(NMAXPI),ZRP(LMAXPI), 3 DALP(NMAXPI,LMAXPI),DALPB(NMAXPI,LMAXPI), 4 DALPC(NMAXPI,LMAXPI),DALPD(NMAXPI,LMAXPI), 5 DRATP(NMAXPI,LMAXPI),DRATPB(NMAXPI,LMAXPI), 6 DRATPC(NMAXPI,LMAXPI),DRATPD(NMAXPI,LMAXPI) real*8 EE(MMAXEI),RE(MMAXEI),REB(MMAXEI),REC(MMAXEI), 1 RED(MMAXEI),YE(MMAXEI),YEB(MMAXEI),YEC(MMAXEI),YED(MMAXEI), 2 TE(NMAXEI),AR(NMAXEI),ARB(NMAXEI),ARC(NMAXEI),ARD(NMAXEI), 3 RS(NMAXEI),RSB(NMAXEI),RSC(NMAXEI),RSD(NMAXEI),BS(LMAXSI), 4 ZRE(LMAXEI),ZS(LMAXTI),ZB(LMAXBI),DALE(NMAXEI,LMAXSI), 5 DALEB(NMAXEI,LMAXSI),DALEC(NMAXEI,LMAXSI),DALED(NMAXEI,LMAXSI), 6 DALB(NMAXEI,LMAXTI),DALBB(NMAXEI,LMAXTI),DALBC(NMAXEI,LMAXTI), 7 DALBD(NMAXEI,LMAXTI),DRATE(NMAXEI,LMAXEI,2), 8 DRATEB(NMAXEI,LMAXEI,2),DRATEC(NMAXEI,LMAXEI,2), 9 DRATED(NMAXEI,LMAXEI,2),DRATB(NMAXEI,LMAXBI,2), X DRATBB(NMAXEI,LMAXBI,2),DRATBC(NMAXEI,LMAXBI,2), 1 DRATBD(NMAXEI,LMAXBI,2) real*8 ZM(IMAXI),Z(IMAXI),ZMM(IMAXI),ZL(IMAXI),TPL(NPTSPI), 1 TPP(NPTSPI),TEL(NPTSEI),ENEWT(NPTSPI),TEE(NPTSEI),RINE(NPTSEI), 2 RINS(NPTSEI),ARES(NPTSEI),YLDE(NPTSEI),DIN(LMAXPI), 3 DINB(LMAXPI),DINC(LMAXPI),DIND(LMAXPI),DRIN(LMAXPI), 4 GP(NPTSPI,IMAXI),GE(NPTSEI,IMAXI,2),GB(NPTSEI,IMAXI,2), 5 EPS(JMAXI),S(JMAXI),SOL(NPTSPI),SPG(NPTSPI),SEG(NPTSEI), 6 G(NPTSI),DOSOL(IMAXI,2),DOSP(IMAXI,2),DOSE(IMAXI,2,2), 7 DOSB(IMAXI,2,2) c C missing declaration real*8 ZMIN,RADCON,ENMU,EMINU,EMAXU,DEP,EMINUL,DELP integer*4 NBEGE,INATT,INEWT,IMIX integer*4 MMAXP,KMAXP,NMAXP,LMAXP,IMIXP integer*4 MMAXE,NMAXE,LMAXS,LMAXE,LMAXT,LMAXB,IMIXE integer*4 M,K,L,N,NLENE,I,J integer*4 NFSTSB,NLSTSB,NLENSB integer*4 NFSTPB,NLSTPB,NLENPB real*8 ANS,RINP,ANSR,BENMU,ZRIN,EMINEL,DEE,DELE integer*4 NE,NP,ISOL,NLENS integer*4 NFSTS,NLSTS,ITRP,NLENP,NFSTP,NLSTP,ILEC,NFSTE,NLSTE real*8 ZRINL,DELTAS,DELTAP,DELTAE,ENEUT,EAV c c declaration added to cleanup the code and avoid memory fault real*8 dosSphIN(IMAXI),dosSpOut(IMAXI) cccc INPUTS: INTEGER*4 IDET, INUC, IMAX, IUNT REAL*8 Zin(IMAXI) REAL*8 EMINS,EMAXS,EMINP,EMAXP,EMINE,EMAXE INTEGER*4 NPTSP,NPTSE INTEGER*4 JSMAX,JPMAX,JEMAX REAL*8 EUNIT,DURATN REAL*8 ESin(JMAXI),SFLUXin(JMAXI) REAL*8 EPin(JMAXI),PFLUXin(JMAXI) REAL*8 EEin(JMAXI),EFLUXin(JMAXI) REAL*8 SolDose(IMAXI,3), ProtDose(IMAXI,3) REAL*8 ElecDose(IMAXI,3), BremDose(IMAXI,3), TotDose(IMAXI,3) ccccccEND OF INPUTS cccccccccccccccccccccccc COMMON /sd2pro/MMAXP,KMAXP,NMAXP,LMAXP,IMIXP COMMON /sd2proTab/EP,RP,TEPN,FEPN,TP,ZRP, &DALP,DRATP COMMON /sd2elbr/MMAXE,NMAXE,LMAXS,LMAXE, &LMAXT,LMAXB,IMIXE COMMON /sd2elbrTab/EE,RE,YE,TE,AR,RS,BS,ZRE &,ZS,ZB,DALE,DALB,DRATE,DRATB DATA DET/'Aluminum','Graphite','Silicon','Air','Bone','CaF2', 1 'GaAs','LiF','SiO2','Tissue','H2O'/ DATA ZMIN/1.0E-06/,RADCON/1.6021892E-08/,NBEGE/1/,ENMU/0.03/ DATA VERSION/'2.10'/ c c add init var to 0 NLSTS = 0 NLSTP = 0 NLSTE = 0 NFSTS = 0 NFSTP = 0 NFSTE = 0 c end add var c C CALL LOGO (VERSION) c PRINT 10 c 10 FORMAT (' Enter input filename: ') c READ 20,FILENM c 20 FORMAT (A) c OPEN (UNIT=9,FILE=FILENM) c READ (9,20) PRTFIL c OPEN (UNIT=10,FILE=PRTFIL) c READ (9,20) ARRFIL c OPEN (UNIT=12,FILE=ARRFIL) c WRITE (10,30) VERSION c WRITE (12,30) VERSION c 30 FORMAT (' OUTPUT FROM SHIELDOSE-2, VERSION ',A) c WRITE (10,32) FILENM c WRITE (12,32) FILENM c 32 FORMAT (' Input filename: ',A) c WRITE (10,34) PRTFIL c WRITE (12,34) PRTFIL c 34 FORMAT (' Print-out filename: ',A) c WRITE (10,36) ARRFIL c WRITE (12,36) ARRFIL c 36 FORMAT (' Array output filename: ',A) c WRITE (10,370) c 370 FORMAT (/' IDET INUC IMAX IUNT') c READ (9,*) IDET,INUC,IMAX,IUNT c WRITE (10,380) IDET,INUC,IMAX,IUNT c 380 FORMAT (12I6) c WRITE (12,380) IDET,IMAX,INUC INATT=2 IF (INUC.EQ.1) INATT=1 INEWT=0 IF (INUC.EQ.3) INEWT=1 ! PRINT *,' Reading database and preparing base arrays.............' ! PRINT *,' Protons.............................................' ! OPEN (UNIT=11,FILE='PROTBAS2.DAT') ! READ (11,20) TAG ! READ (11,*) MMAXP,KMAXP,NMAXP,LMAXP,IMIX ! READ (11,*) (EP(M),M=1,MMAXP) ! READ (11,*) (RP(M),M=1,MMAXP) ! READ (11,*) (TEPN(K),K=1,KMAXP) ! READ (11,*) (FEPN(K),K=1,KMAXP) ! READ (11,*) (TP(N),N=1,NMAXP) ! READ (11,*) (ZRP(L),L=1,LMAXP) ! DO 50 N=1,NMAXP ! DO 38 I=1,2 ! READ (11,*) (DUM(L),L=1,LMAXP) ! IF (I.NE.INATT) GO TO 38 ! DO 390 L=1,LMAXP ! 390 DALP(N,L)=DUM(L) ! DALP(N,LMAXP)=(DALP(N,LMAXP-1)/DALP(N,LMAXP-2))*DALP(N,LMAXP-1) ! 38 CONTINUE ! DO 40 I=1,IMIX ! READ (11,*) (DUM(L),L=1,LMAXP) ! IF (I.NE.IDET) GO TO 40 ! DO 39 L=1,LMAXP ! 39 DRATP(N,L)=DUM(L) ! 40 CONTINUE ! 50 CONTINUE ! CLOSE (11) call protbas2(INATT,IDET,tag) IMIX=IMIXP DO 60 M=1,MMAXP EP(M)=LOG(EP(M)) 60 RP(M)=LOG(RP(M)) DO 65 K=1,KMAXP 65 TEPN(K)=LOG(TEPN(K)) DO 70 N=1,NMAXP 70 TP(N)=LOG(TP(N)) CALL SCOF (MMAXP,EP,RP,RPB,RPC,RPD) CALL SCOF (KMAXP,TEPN,FEPN,FEPNB,FEPNC,FEPND) DO 75 L=1,LMAXP DO 72 N=1,NMAXP 72 DALP(N,L)=LOG(DALP(N,L)) CALL SCOF (NMAXP,TP,DALP(1,L),DALPB(1,L),DALPC(1,L),DALPD(1,L)) 75 CONTINUE ! PRINT *,' Electrons and bremsstrahlung........................' ! OPEN (UNIT=11,FILE='ELBRBAS2.DAT') ! READ (11,20) TAG ! READ (11,*) MMAXE,NMAXE,LMAXS,LMAXE,LMAXT,LMAXB,IMIX call elbrbas2(IDET,tag) IMIX=IMIXE NLENE=NMAXE-NBEGE+1 ! READ (11,*) (EE(M),M=1,MMAXE) ! READ (11,*) (RE(M),M=1,MMAXE) ! READ (11,*) (YE(M),M=1,MMAXE) ! READ (11,*) (TE(N),N=1,NMAXE) ! READ (11,*) (AR(N),N=1,NMAXE) ! READ (11,*) (RS(N),N=1,NMAXE) ! READ (11,*) (BS(L),L=1,LMAXS) ! BS(LMAXS+1)=2.0 ! READ (11,*) (ZRE(L),L=1,LMAXE) ! READ (11,*) (ZS(L),L=1,LMAXT) ! READ (11,*) (ZB(L),L=1,LMAXB) ! DO 100 N=1,NMAXE ! READ (11,*) (DALE(N,L),L=1,LMAXS) ! DALE(N,LMAXS+1)=1.0E-07 ! READ (11,*) (DALB(N,L),L=1,LMAXT) ! DO 90 I=1,IMIX ! DO 80 M=1,2 ! READ (11,*) (DUM(L),L=1,LMAXE) ! IF (I.NE.IDET) GO TO 77 ! DO 76 L=1,LMAXE ! 76 DRATE(N,L,M)=DUM(L) ! 77 READ (11,*) (DUM(L),L=1,LMAXB) ! IF (I.NE.IDET) GO TO 80 ! DO 78 L=1,LMAXB ! 78 DRATB(N,L,M)=DUM(L) ! 80 CONTINUE ! 90 CONTINUE ! 100 CONTINUE LMAXS=LMAXS+1 ! CLOSE (11) DO 110 M=1,MMAXE EE(M)=LOG(EE(M)) RE(M)=LOG(RE(M)) 110 YE(M)=LOG(YE(M)) DO 120 N=1,NMAXE TE(N)=LOG(TE(N)) AR(N)=LOG(AR(N)) 120 RS(N)=LOG(RS(N)) DO 130 L=1,LMAXB 130 ZB(L)=LOG(ZB(L)) CALL SCOF (MMAXE,EE,RE,REB,REC,RED) CALL SCOF (MMAXE,EE,YE,YEB,YEC,YED) CALL SCOF (NMAXE,TE,AR,ARB,ARC,ARD) CALL SCOF (NMAXE,TE,RS,RSB,RSC,RSD) DO 150 L=1,LMAXS DO 140 N=NBEGE,NMAXE 140 DALE(N,L)=LOG(DALE(N,L)) CALL LCOF (NLENE,TE(NBEGE),DALE(NBEGE,L),DALEB(NBEGE,L), 1 DALEC(NBEGE,L),DALED(NBEGE,L)) C CALL SCOF (NLENE,TE(NBEGE),DALE(NBEGE,L),DALEB(NBEGE,L), C 1 DALEC(NBEGE,L),DALED(NBEGE,L)) 150 CONTINUE DO 170 L=1,LMAXT ZS(L)=LOG(ZS(L)) DO 160 N=NBEGE,NMAXE 160 DALB(N,L)=LOG(DALB(N,L)) CALL LCOF (NLENE,TE(NBEGE),DALB(NBEGE,L),DALBB(NBEGE,L), 1 DALBC(NBEGE,L),DALBD(NBEGE,L)) C CALL SCOF (NLENE,TE(NBEGE),DALB(NBEGE,L),DALBB(NBEGE,L), C 1 DALBC(NBEGE,L),DALBD(NBEGE,L)) 170 CONTINUE ! PRINT *,' Preparing base arrays for selected detector material...' DO 220 L=1,LMAXP CALL SCOF (NMAXP,TP,DRATP(1,L),DRATPB(1,L),DRATPC(1,L), 1 DRATPD(1,L)) 220 CONTINUE DO 240 M=1,2 DO 230 L=1,LMAXE CALL LCOF (NLENE,TE(NBEGE),DRATE(NBEGE,L,M), 1 DRATEB(NBEGE,L,M),DRATEC(NBEGE,L,M),DRATED(NBEGE,L,M)) C CALL SCOF (NLENE,TE(NBEGE),DRATE(NBEGE,L,M), C 1 DRATEB(NBEGE,L,M),DRATEC(NBEGE,L,M),DRATED(NBEGE,L,M)) 230 CONTINUE 240 CONTINUE DO 260 M=1,2 DO 250 L=1,LMAXB CALL LCOF (NLENE,TE(NBEGE),DRATB(NBEGE,L,M), 1 DRATBB(NBEGE,L,M),DRATBC(NBEGE,L,M),DRATBD(NBEGE,L,M)) C CALL SCOF (NLENE,TE(NBEGE),DRATB(NBEGE,L,M), C 1 DRATBB(NBEGE,L,M),DRATBC(NBEGE,L,M),DRATBD(NBEGE,L,M)) 250 CONTINUE 260 CONTINUE GO TO (440,470,500), IUNT 440 continue c 440 WRITE (10,450) c 450 FORMAT (/' SHIELD DEPTH (mils)') c READ (9,*) (ZM(I),I=1,IMAX) c WRITE (10,455) (ZM(I),I=1,IMAX) c 455 FORMAT (1P6E12.5) do I = 1,IMAX ZM(I) = Zin(I) enddo DO 460 I=1,IMAX IF (ZM(I).LE.ZMIN/ZCON) ZM(I)=ZMIN/ZCON Z(I)=ZCON*ZM(I) 460 ZMM(I)=Z(I)*ZMCON GO TO 530 470 continue c 470 WRITE (10,480) c 480 FORMAT (/' SHIELD DEPTH (g/cm2)') c READ (9,*) (Z(I),I=1,IMAX) c WRITE (10,455) (Z(I),I=1,IMAX) DO I=1,IMAX Z(I) = Zin(I) ENDDO DO 490 I=1,IMAX IF (Z(I).LE.ZMIN) Z(I)=ZMIN ZM(I)=Z(I)/ZCON 490 ZMM(I)=Z(I)*ZMCON GO TO 530 500 continue c 500 WRITE (10,510) c 510 FORMAT (/' SHIELD DEPTH (mm)') c READ (9,*) (ZMM(I),I=1,IMAX) c WRITE (10,455) (ZMM(I),I=1,IMAX) DO I=1,IMAX ZMM(I) = Zin(I) ENDDO DO 520 I=1,IMAX IF (ZMM(I).LE.ZMIN*ZMCON) ZMM(I)=ZMIN*ZMCON Z(I)=ZMM(I)/ZMCON 520 ZM(I)=Z(I)/ZCON 530 DO 540 I=1,IMAX 540 ZL(I)=LOG(Z(I)) c WRITE (12,1435) (Z(I),I=1,IMAX) c WRITE (10,550) c 550 FORMAT(/' EMINS EMAXS EMINP EMAXP NPTSP EMINE c 1 EMAXE NPTSE') c READ (9,*) EMINS,EMAXS,EMINP,EMAXP,NPTSP,EMINE,EMAXE,NPTSE c WRITE (10,560) EMINS,EMAXS,EMINP,EMAXP,NPTSP,EMINE,EMAXE,NPTSE c 560 FORMAT (4F10.3,I6,2F10.3,I6) EMINU=MIN(EMINP,EMINS) EMAXU=MAX(EMAXP,EMAXS) DEP=LOG(EMAXU/EMINU)/FLOAT(NPTSP-1) EMINUL=LOG(EMINU) DELP=DEP/3.0 CALL EINDEX (EMINU,DEP,NPTSP,EMINS,EMAXS,NFSTSB,NLSTSB,NLENSB) CALL EINDEX (EMINU,DEP,NPTSP,EMINP,EMAXP,NFSTPB,NLSTPB,NLENPB) DO 570 NP=1,NPTSP TPL(NP)=EMINUL+FLOAT(NP-1)*DEP TPP(NP)=EXP(TPL(NP)) 570 CONTINUE c WRITE (10,580) TPP(NFSTSB),TPP(NLSTSB),TPP(NFSTPB),TPP(NLSTPB), c 1 NPTSP,EMINE,EMAXE,NPTSE c 580 FORMAT (4F10.3,I6,2F10.3,I6,' ADJUSTED VALUES') c WRITE (12,560) TPP(NFSTSB),TPP(NLSTSB),TPP(NFSTPB),TPP(NLSTPB), c 1 NPTSP,EMINE,EMAXE,NPTSE ! PRINT *,' Preparing mesh arrays to be integrated over spectra....' ! PRINT *,' Protons.............................................' DO 660 NP=1,NPTSP CALL BSPOL (TPL(NP),MMAXP,EP,RP,RPB,RPC,RPD,ANS) RINP=EXP(ANS) DO 610 L=1,LMAXP IF (TPL(NP).LT.TP(NMAXP)) GO TO 590 ANS=DALP(NMAXP,L) ANSR=DRATP(NMAXP,L) GO TO 605 590 IF (TPL(NP).GT.TP(1)) GO TO 600 ANS=DALP(1,L) ANSR=DRATP(1,L) GO TO 605 600 CALL BSPOL (TPL(NP),NMAXP,TP,DALP(1,L),DALPB(1,L), 1 DALPC(1,L),DALPD(1,L),ANS) ANSR=1.0 IF (IDET.EQ.1) GO TO 605 CALL BSPOL (TPL(NP),NMAXP,TP,DRATP(1,L),DRATPB(1,L), 1 DRATPC(1,L),DRATPD(1,L),ANSR) 605 DIN(L)=ANS+LOG(ANSR) 610 CONTINUE ENEWT(NP)=0.0 BENMU=0.0 IF (INATT.EQ.1) GO TO 620 IF (TPL(NP).LE.TEPN(1)) GO TO 615 CALL BSPOL (TPL(NP),KMAXP,TEPN,FEPN,FEPNB,FEPNC,FEPND,ANS) ENEWT(NP)=TPP(NP)*ANS 615 BENMU=ENEWT(NP)*ENMU 620 CALL SCOF (LMAXP,ZRP,DIN,DINB,DINC,DIND) DO 650 I=1,IMAX ZRIN=Z(I)/RINP IF (ZRIN.LT.ZRP(LMAXP)) GO TO 640 GP(NP,I)=0.0 GO TO 645 640 CALL BSPOL (ZRIN,LMAXP,ZRP,DIN,DINB,DINC,DIND,ANS) ANS=EXP(ANS) GP(NP,I)=TPP(NP)*ANS/RINP 645 IF (INEWT.EQ.1.AND.TPL(NP).GT.TEPN(1)) GP(NP,I)=GP(NP,I)+ 1 BENMU*EXP(-ENMU*Z(I)) 650 CONTINUE 660 CONTINUE ! PRINT *,' Electrons and bremsstrahlung try...................' EMINEL=LOG(EMINE) DEE=(LOG(EMAXE)-EMINEL)/FLOAT(NPTSE-1) DELE=DEE/3.0 DO 670 NE=1,NPTSE TEL(NE)=EMINEL+FLOAT(NE-1)*DEE TEE(NE)=EXP(TEL(NE)) CALL BSPOL (TEL(NE),MMAXE,EE,RE,REB,REC,RED,ANS) RINE(NE)=EXP(ANS) CALL BSPOL (TEL(NE),NMAXE,TE,RS,RSB,RSC,RSD,ANS) RINS(NE)=RINE(NE)*EXP(ANS) CALL BSPOL (TEL(NE),NMAXE,TE,AR,ARB,ARC,ARD,ANS) ARES(NE)=EXP(ANS) CALL BSPOL (TEL(NE),MMAXE,EE,YE,YEB,YEC,YED,ANS) 670 YLDE(NE)=EXP(ANS) DO 820 M=1,2 DO 815 NE=1,NPTSE DO 700 L=1,LMAXS IF (TEL(NE).LT.TE(NMAXE)) GO TO 680 DIN(L)=DALE(NMAXE,L) GO TO 700 680 IF (TEL(NE).GT.TE(NBEGE)) GO TO 690 DIN(L)=DALE(NBEGE,L) GO TO 700 690 CALL BSPOL (TEL(NE),NLENE,TE(NBEGE),DALE(NBEGE,L),DALEB(NBEGE,L), 1 DALEC(NBEGE,L),DALED(NBEGE,L),DIN(L)) 700 CONTINUE DO 715 L=1,LMAXE DRIN(L)=1.0 IF (IDET.EQ.1.AND.M.EQ.1) GO TO 715 IF (TEL(NE).LT.TE(NMAXE)) GO TO 710 DRIN(L)=DRATE(NMAXE,L,M) GO TO 715 710 IF (TEL(NE).GT.TE(NBEGE)) GO TO 712 DRIN(L)=DRATE(NBEGE,L,M) GO TO 715 712 CALL BSPOL (TEL(NE),NLENE,TE(NBEGE),DRATE(NBEGE,L,M), 1 DRATEB(NBEGE,L,M),DRATEC(NBEGE,L,M),DRATED(NBEGE,L,M),DRIN(L)) IF (DRIN(L).LT.0.0) DRIN(L)=0.0 715 CONTINUE C CALL LCOF (LMAXS,BS,DIN,DINB,DINC,DIND) CALL SCOF (LMAXS,BS,DIN,DINB,DINC,DIND) DO 740 I=1,IMAX ZRIN=Z(I)/RINS(NE) IF (ZRIN.LT.BS(LMAXS)) GO TO 730 720 GE(NE,I,M)=0.0 GO TO 740 730 CALL BSPOL (ZRIN,LMAXS,BS,DIN,DINB,DINC,DIND,ANS) ANS=EXP(ANS) GE(NE,I,M)=TEE(NE)*ANS*ARES(NE)/RINS(NE) 740 CONTINUE C CALL LCOF (LMAXE,ZRE,DRIN,DINB,DINC,DIND) CALL SCOF (LMAXE,ZRE,DRIN,DINB,DINC,DIND) DO 745 I=1,IMAX ZRIN=Z(I)/RINE(NE) IF (ZRIN.LT.ZRE(LMAXE)) GO TO 742 GE(NE,I,M)=GE(NE,I,M)*DRIN(LMAXE) GO TO 745 742 CALL BSPOL (ZRIN,LMAXE,ZRE,DRIN,DINB,DINC,DIND,ANSR) IF (ANSR.LT.0.0) ANSR=0.0 GE(NE,I,M)=GE(NE,I,M)*ANSR 745 CONTINUE DO 780 L=1,LMAXT IF (TEL(NE).LT.TE(NMAXE)) GO TO 760 DIN(L)=DALB(NMAXE,L) GO TO 780 760 IF (TEL(NE).GT.TE(NBEGE)) GO TO 770 DIN(L)=DALB(NBEGE,L) GO TO 780 770 CALL BSPOL (TEL(NE),NLENE,TE(NBEGE),DALB(NBEGE,L),DALBB(NBEGE,L), 1 DALBC(NBEGE,L),DALBD(NBEGE,L),DIN(L)) 780 CONTINUE DO 795 L=1,LMAXB DRIN(L)=1.0 IF (IDET.EQ.1.AND.M.EQ.1) GO TO 795 IF (TEL(NE).LT.TE(NMAXE)) GO TO 790 DRIN(L)=DRATB(NMAXE,L,M) GO TO 795 790 IF (TEL(NE).GT.TE(NBEGE)) GO TO 792 DRIN(L)=DRATB(NBEGE,L,M) GO TO 795 792 CALL BSPOL (TEL(NE),NLENE,TE(NBEGE),DRATB(NBEGE,L,M), 1 DRATBB(NBEGE,L,M),DRATBC(NBEGE,L,M),DRATBD(NBEGE,L,M),DRIN(L)) IF (DRIN(L).LT.0.0) DRIN(L)=0.0 795 CONTINUE CALL LCOF (LMAXT,ZS,DIN,DINB,DINC,DIND) C CALL SCOF (LMAXT,ZS,DIN,DINB,DINC,DIND) DO 800 I=1,IMAX ZRINL=LOG(Z(I)/RINE(NE)) CALL BSPOL (ZRINL,LMAXT,ZS,DIN,DINB,DINC,DIND,ANS) ANS=EXP(ANS) GB(NE,I,M)=TEE(NE)*ANS*YLDE(NE)/RINE(NE) 800 CONTINUE CALL LCOF (LMAXB,ZB,DRIN,DINB,DINC,DIND) C CALL SCOF (LMAXB,ZB,DRIN,DINB,DINC,DIND) DO 812 I=1,IMAX IF (ZL(I).LT.ZB(LMAXB)) GO TO 810 GB(NE,I,M)=GB(NE,I,M)*DRIN(LMAXB) GO TO 812 810 CALL BSPOL (ZL(I),LMAXB,ZB,DRIN,DINB,DINC,DIND,ANSR) IF (ANSR.LT.0.0) ANSR=0.0 GB(NE,I,M)=GB(NE,I,M)*ANSR 812 CONTINUE 815 CONTINUE 820 CONTINUE ! PRINT *,' Performing calculations for input spectra..............' 830 continue c 830 WRITE (10,840) c 840 FORMAT (/) c WRITE (10,850) c 850 FORMAT (' ') c READ (9,20,END=1440) TAG c PRINT 860, TAG c 860 FORMAT (4X,A) c WRITE (10,20) TAG c WRITE (12,20) TAG c WRITE (10,870) c 870 FORMAT (/' JSMAX JPMAX JEMAX EUNIT DURATN') c READ (9,*) JSMAX,JPMAX,JEMAX,EUNIT,DURATN c WRITE (10,880) JSMAX,JPMAX,JEMAX,EUNIT,DURATN c WRITE (12,880) JSMAX,JPMAX,JEMAX,EUNIT,DURATN c 880 FORMAT (3I6,1P2E12.5) IF (DURATN.LE.0.0) DURATN=1.0 DELTAS=RADCON*DELP/4.0 DELTAP=DURATN*RADCON*DELP/4.0 DELTAE=DURATN*RADCON*DELE/4.0 IF (EUNIT.LE.0.0) EUNIT=1.0 ISOL=2 IF (JSMAX.LT.3) GO TO 900 ISOL=1 c WRITE (10,885) c 885 FORMAT (//' E(MeV)') c READ (9,*) (EPS(J),J=1,JSMAX) DO J=1,JSMAX EPS(J) = ESin(J) ENDDO c WRITE (10,905) (EPS(J),J=1,JSMAX) c WRITE (12,905) (EPS(J),J=1,JSMAX) c WRITE (10,890) c 890 FORMAT (/' SOLAR PROTON SPECTRUM (/energy/cm2)') c READ (9,*) (S(J),J=1,JSMAX) DO J=1,JSMAX S(J) = SFLUXin(J) ENDDO c WRITE (10,905) (S(J),J=1,JSMAX) c WRITE (12,905) (S(J),J=1,JSMAX) NLENS=NLENSB NFSTS=NFSTSB NLSTS=NLSTSB CALL SPECTR (JSMAX,EPS,S,EUNIT,EMINU,DEP,NPTSP,NFSTS,NLSTS,NLENS, 1 TPP,TPL,SOL) c WRITE (10,891) TPP(NFSTS),TPP(NLSTS),NLENS c 891 FORMAT (/' SPECTRUM INTEGRATED FROM',1PE11.4,' TO',1PE11.4, c 1 ' MeV, USING',I5,' POINTS') DO 892 NP=NFSTS,NLSTS 892 G(NP)=SOL(NP)*ENEWT(NP) CALL INTEG (DELP,G(NFSTS),NLENS,ENEUT) DO 894 NP=NFSTS,NLSTS 894 G(NP)=SOL(NP)*TPP(NP) CALL INTEG (DELP,G(NFSTS),NLENS,EAV) ENEUT=ENEUT/EAV c WRITE (10,896) ENEUT c 896 FORMAT (/' ASSUMED FRACTION OF BEAM ENERGY INTO NEUTRON ENERGY =', c 1 1PE12.5) 900 ITRP=2 IF (JPMAX.LT.3) GO TO 920 ITRP=1 c WRITE (10,885) c READ (9,*) (EPS(J),J=1,JPMAX) DO J=1,JPMAX EPS(J) = EPin(J) ENDDO c WRITE (10,905) (EPS(J),J=1,JPMAX) c WRITE (12,905) (EPS(J),J=1,JPMAX) c 905 FORMAT (1P10E12.4) c WRITE (10,910) c 910 FORMAT (/' TRAPPED PROTON SPECTRUM (/energy/cm2/time)') c READ (9,*) (S(J),J=1,JPMAX) DO J=1,JPMAX S(J) = PFLUXin(J) ENDDO c WRITE (10,905) (S(J),J=1,JPMAX) c WRITE (12,905) (S(J),J=1,JPMAX) NLENP=NLENPB NFSTP=NFSTPB NLSTP=NLSTPB CALL SPECTR (JPMAX,EPS,S,EUNIT,EMINU,DEP,NPTSP,NFSTP,NLSTP,NLENP, 1 TPP,TPL,SPG) c WRITE (10,891) TPP(NFSTP),TPP(NLSTP),NLENP DO 912 NP=NFSTP,NLSTP 912 G(NP)=SPG(NP)*ENEWT(NP) CALL INTEG (DELP,G(NFSTP),NLENP,ENEUT) DO 914 NP=NFSTP,NLSTP 914 G(NP)=SPG(NP)*TPP(NP) CALL INTEG (DELP,G(NFSTP),NLENP,EAV) ENEUT=ENEUT/EAV c WRITE (10,896) ENEUT 920 ILEC=2 IF (JEMAX.LT.3) GO TO 940 ILEC=1 c WRITE (10,885) c READ (9,*) (EPS(J),J=1,JEMAX) DO J=1,JEMAX EPS(J) = EEin(J) ENDDO c WRITE (10,905) (EPS(J),J=1,JEMAX) c WRITE (12,905) (EPS(J),J=1,JEMAX) c WRITE (10,930) c 930 FORMAT (/' ELECTRON SPECTRUM (/energy/cm2/time)') c READ (9,*) (S(J),J=1,JEMAX) DO J=1,JEMAX S(J) = EFLUXin(J) ENDDO c WRITE (10,905) (S(J),J=1,JEMAX) c WRITE (12,905) (S(J),J=1,JEMAX) NLENE=NPTSE NFSTE=1 NLSTE=NPTSE CALL SPECTR (JEMAX,EPS,S,EUNIT,EMINE,DEE,NPTSE,NFSTE,NLSTE,NLENE, 1 TEE,TEL,SEG) c c WRITE (10,891) TEE(NFSTE),TEE(NLSTE),NLENE 940 GO TO (980,950), ISOL ! 950 DO 960 NP=NFSTS,NLSTS ! 960 SOL(NP)=0.0 950 continue DO 970 J=1,2 DO 970 I=1,IMAX 970 DOSOL(I,J)=0.0 GO TO 1010 980 DO 1000 I=1,IMAX DO 990 NP=NFSTS,NLSTS 990 G(NP)=SOL(NP)*GP(NP,I) CALL INTEG (DELTAS,G(NFSTS),NLENS,DOSOL(I,1)) 1000 CONTINUE do I=1,IMAX dosSphIN(I)=DOSOL(I,1) enddo CALL SPHERE (ZL,dosSphIN,IMAX,dosSpOut) do I=1,IMAX DOSOL(I,2)=dosSpOut(I) enddo 1010 GO TO (1050,1020), ITRP c 1020 DO 1030 NP=NFSTP,NLSTP c 1030 SPG(NP)=0.0 1020 continue DO 1040 J=1,2 DO 1040 I=1,IMAX 1040 DOSP(I,J)=0.0 GO TO 1080 1050 DO 1070 I=1,IMAX DO 1060 NP=NFSTP,NLSTP 1060 G(NP)=SPG(NP)*GP(NP,I) CALL INTEG (DELTAP,G(NFSTP),NLENP,DOSP(I,1)) 1070 CONTINUE do I=1,IMAX dosSphIN(I)=DOSP(I,1) enddo CALL SPHERE (ZL,dosSphIN,IMAX,dosSpOut) do I=1,IMAX DOSP(I,2)=dosSpOut(I) enddo 1080 GO TO (1110,1090), ILEC 1090 DO 1100 J=1,2 DO 1100 M=1,2 DO 1100 I=1,IMAX DOSE(I,M,J)=0.0 1100 DOSB(I,M,J)=0.0 GO TO 1160 1110 DO 1150 M=1,2 DO 1130 I=1,IMAX DO 1120 NE=NFSTE,NLSTE G(NE)=SEG(NE)*GE(NE,I,M) 1120 SPG(NE)=SEG(NE)*GB(NE,I,M) CALL INTEG (DELTAE,G(NFSTE),NLENE,DOSE(I,M,1)) CALL INTEG (DELTAE,SPG(NFSTE),NLENE,DOSB(I,M,1)) 1130 CONTINUE GO TO (1140,1150), M 1140 continue do I=1,IMAX dosSphIN(I)=DOSE(I,M,1) enddo CALL SPHERE (ZL,dosSphIN,IMAX,dosSpOut) do I=1,IMAX DOSE(I,M,2)=dosSpOut(I) enddo do I=1,IMAX dosSphIN(I)=DOSB(I,M,1) enddo CALL SPHERE (ZL,dosSphIN,IMAX,dosSpOut) do I=1,IMAX DOSB(I,M,2)=dosSpOut(I) enddo 1150 CONTINUE 1160 continue ccccccccccccccccccccc my code to store outputs c K=1, J=1,M=1, DOSE IN SEMI-INFINITE ALUMINUM MEDIUM c K=2, J=1,M=2, DOSE AT TRANSMISSION SURFACE OF FINITE ALUMINUM SLAB SHIELDS c K=3, J=2,M=1, 1/2 DOSE AT CENTER OF ALUMINUM SPHERES c DOSOL(I,J), DOSP(I,J), DOSE(I,J,M), DOSB(I,J,M) ! write(6,*)'Format output ...' K=1 J=1 M=1 DO I=1,IMAX ProtDose(I,K) = DOSP(I,J) SolDose(I,K) = DOSOL(I,J) ElecDose(I,K) = DOSE(I,M,J) BremDose(I,K) = DOSB(I,M,J) TotDose(I,K)=DOSP(I,J)+DOSOL(I,J)+DOSE(I,M,J)+DOSB(I,M,J) ENDDO K=2 J=1 M=2 DO I=1,IMAX ProtDose(I,K) = DOSP(I,J) SolDose(I,K) = DOSOL(I,J) ElecDose(I,K) = DOSE(I,M,J) BremDose(I,K) = DOSB(I,M,J) TotDose(I,K)=DOSP(I,J)+DOSOL(I,J)+DOSE(I,M,J)+DOSB(I,M,J) ENDDO K=3 J=2 M=1 DO I=1,IMAX ProtDose(I,K) = DOSP(I,J) SolDose(I,K) = DOSOL(I,J); ElecDose(I,K) = DOSE(I,M,J) BremDose(I,K) = DOSB(I,M,J) TotDose(I,K)=DOSP(I,J)+DOSOL(I,J)+DOSE(I,M,J)+DOSB(I,M,J) ENDDO c 1160 return ccccccccccccccccccccc all the rest of this code is just write statements c 1160 J=1 c DO 1340 M=2,1,-1 c GO TO (1190,1170), M c 1170 WRITE (10,1180) c 1180 FORMAT(//' DOSE AT TRANSMISSION SURFACE OF FINITE ALUMINUM SLAB SH c 1IELDS') c GO TO 1210 c 1190 WRITE (10,1200) c 1200 FORMAT(//' DOSE IN SEMI-INFINITE ALUMINUM MEDIUM') c 1210 WRITE (10,1230) DET(IDET) c 1230 FORMAT (/' rads ',A) c IF (INATT.EQ.1) WRITE (10,1240) c 1240 FORMAT (/' Proton results without nuclear attenuation') c IF (INATT.EQ.2) WRITE (10,1250) c 1250 FORMAT (/' Proton results with approximate treatment of nuclear at c 1tenuation') c IF (INATT.EQ.2.AND.INEWT.EQ.0) WRITE (10,1260) c 1260 FORMAT ( ' neglecting transport of energy by neutrons') c IF (INATT.EQ.2.AND.INEWT.EQ.1) WRITE (10,1270) c 1270 FORMAT ( ' and crude exponential transport of energy by neutron c 1s') c WRITE (10,1310) c 1310 FORMAT(/' Z(mils) Z(mm) Z(g/cm2) ELECTRON BREMS c 1 EL+BR TRP PROT SOL PROT EL+BR+TRP TOTAL') c WRITE (10,850) c DO 1330 I=1,IMAX c DOSEB=DOSE(I,M,J)+DOSB(I,M,J) c DOSEBP=DOSEB+DOSP(I,J) c DOST=DOSEBP+DOSOL(I,J) c WRITE (10,1320) ZM(I),ZMM(I),Z(I),DOSE(I,M,J),DOSB(I,M,J),DOSEB, c 1 DOSP(I,J),DOSOL(I,J),DOSEBP,DOST c 1320 FORMAT (1P10E11.3) c IF (FLOAT(I/10).EQ.0.1*FLOAT(I)) WRITE (10,850) c 1330 CONTINUE c 1340 CONTINUE c J=2 c M=1 c WRITE (10,1350) c 1350 FORMAT (//' 1/2 DOSE AT CENTER OF ALUMINUM SPHERES') c WRITE (10,1230) DET(IDET) c IF (INATT.EQ.1) WRITE (10,1240) c IF (INATT.EQ.2) WRITE (10,1250) c IF (INATT.EQ.2.AND.INEWT.EQ.0) WRITE (10,1260) c IF (INATT.EQ.2.AND.INEWT.EQ.1) WRITE (10,1270) c WRITE (10,1310) c WRITE (10,850) c DO 1410 I=1,IMAX c DOSEB=DOSE(I,M,J)+DOSB(I,M,J) c DOSEBP=DOSEB+DOSP(I,J) c DOST=DOSEBP+DOSOL(I,J) c WRITE (10,1320) ZM(I),ZMM(I),Z(I),DOSE(I,M,J),DOSB(I,M,J),DOSEB, c 1 DOSP(I,J),DOSOL(I,J),DOSEBP,DOST c IF (FLOAT(I/10).EQ.0.1*FLOAT(I)) WRITE (10,850) c 1410 CONTINUE c DO 1437 J=1,2 c WRITE (12,1435) (DOSOL(I,J),I=1,IMAX) c WRITE (12,1435) (DOSP(I,J),I=1,IMAX) c 1435 FORMAT (1P10E10.3) c 1437 CONTINUE c DO 1438 M=2,1,-1 c WRITE (12,1435) (DOSE(I,M,1),I=1,IMAX) c WRITE (12,1435) (DOSB(I,M,1),I=1,IMAX) c 1438 CONTINUE c WRITE (12,1435) (DOSE(I,1,2),I=1,IMAX) c WRITE (12,1435) (DOSB(I,1,2),I=1,IMAX) c GO TO 830 c 1440 PRINT 32, FILENM c PRINT 34, PRTFIL c PRINT 36, ARRFIL cC print 1500,CHAR(27) cC 1500 format (' ',A1,'[0m') c STOP END C SUBROUTINE EINDEX (EMINB,DE,NPTS,EMIN,EMAX,NFST,NLST,NLEN), 28 APR 94. SUBROUTINE EINDEX (EMINB,DE,NPTS,EMIN,EMAX,NFST,NLST,NLEN) implicit none C IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer*4 NFST,NLST,NLEN,NPTS real*8 EMINB,EMIN,EMAX,DE NFST=LOG(EMIN/EMINB)/DE+0.5 NFST=NFST+1 IF (NFST.LT.1) NFST=1 NLST=LOG(EMAX/EMINB)/DE+0.5 NLST=NLST+1 IF (NLST.GT.NPTS) NLST=NPTS NLEN=NLST-NFST+1 RETURN END C SUBROUTINE SPECTR (JMAX,EPS,S,EUNIT,EMINB,DEL,NPTS,NFST,NLST,NLEN, C 1 T,TL,SP), 28 APR 94. SUBROUTINE SPECTR (JMAX,EPS,S,EUNIT,EMINB,DEL,NPTS,NFST,NLST,NLEN, 1 T,TL,SP) implicit none c IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer*4 JMAXI,NPTSI PARAMETER (JMAXI=301,NPTSI=1001) integer*4 JMAX,NPTS,NFST,NLST,NLEN,N,J real*8 EPS(JMAXI),S(JMAXI),T(NPTSI),TL(NPTSI),SP(NPTSI), & BCOF(JMAXI),CCOF(JMAXI), 1 DCOF(JMAXI),G(NPTSI) real*8 EUNIT,EMINB,P,ARG real*8 EMC2,ARGLIM,DELTA,DEL,ALPHA,BETA,ANS,SIN,EBAR C ARGLIM = 700.0 (DOUBLE PRECISION), 85.0 (SINGLE PRECISION) DATA EMC2/938.27231/,ARGLIM/85.0/ DELTA=DEL/3.0 IF (EPS(1).GT.0.0) GO TO 20 ALPHA=S(1) BETA=S(2) IF (BETA.LE.0.0) BETA=1.0 BETA=BETA/ALPHA DO 10 N=NFST,NLST SP(N)=0.0 G(N)=0.0 IF (S(3).LE.0.0) GO TO 6 P=SQRT(T(N)*(T(N)+2.0*EMC2)) ARG=P/ALPHA IF (ARG.GT.ARGLIM) GO TO 10 SP(N)=T(N)*BETA*((T(N)+EMC2)/P)*EXP(-ARG) GO TO 8 6 ARG=T(N)/ALPHA IF (ARG.GT.ARGLIM) GO TO 10 SP(N)=T(N)*BETA*EXP(-ARG) 8 G(N)=T(N)*SP(N) 10 CONTINUE GO TO 50 20 CALL EINDEX (EMINB,DEL,NPTS,EPS(1),EPS(JMAX),NFST,NLST,NLEN) DO 30 J=1,JMAX EPS(J)=LOG(EPS(J)) 30 S(J)=LOG(EUNIT*S(J)) CALL SCOF (JMAX,EPS,S,BCOF,CCOF,DCOF) DO 40 N=NFST,NLST CALL BSPOL (TL(N),JMAX,EPS,S,BCOF,CCOF,DCOF,ANS) SP(N)=T(N)*EXP(ANS) 40 G(N)=T(N)*SP(N) 50 CALL INTEG (DELTA,SP(NFST),NLEN,SIN) CALL INTEG (DELTA,G(NFST),NLEN,EBAR) EBAR=EBAR/SIN c WRITE (10,60) c 60 FORMAT (/' INT SPEC EAV(MeV)') c WRITE (10,70) SIN,EBAR c 70 FORMAT (1PE12.4,0PF12.5) RETURN END C SUBROUTINE SPHERE (ZL,DOSE,IMAX,DOSPH), 26 JAN 93. SUBROUTINE SPHERE (ZL,DOSE,IMAX,DOSPH) implicit none C IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer*4 IMAXI,IMAX,I,IMIX,IMIX1 PARAMETER (IMAXI=71) real*8 ZL(IMAX),DOSE(IMAXI),DOSPH(IMAXI),DOSL(IMAXI),BCOF(IMAXI), 1 CCOF(IMAXI),DCOF(IMAXI) DO 10 I=1,IMAX IF (DOSE(I).LE.0.0) GO TO 20 10 DOSL(I)=LOG(DOSE(I)) I=IMAX+1 20 IMIX=I-1 IF (IMIX.LT.3) GO TO 40 CALL SCOF (IMIX,ZL,DOSL,BCOF,CCOF,DCOF) BCOF(IMIX)=BCOF(IMIX-1)+(2.0*CCOF(IMIX-1)+3.0*DCOF(IMIX-1)* 1 (ZL(IMIX)-ZL(IMIX-1)))*(ZL(IMIX)-ZL(IMIX-1)) DO 30 I=1,IMIX 30 DOSPH(I)=DOSE(I)*(1.0-BCOF(I)) 40 IMIX1=IMIX+1 IF (IMIX1.GT.IMAX) RETURN DO 50 I=IMIX1,IMAX 50 DOSPH(I)=0.0 RETURN END SUBROUTINE SCOF(N,X,Y,B,C,D) implicit none C REINSCH ALGORITHM, VIA MJB, 22 FEB 83 C Y(S)=((D(J)*(X-X(J))+C(J))*(X-X(J))+B(J))*(X-X(J))+Y(J) C FOR X BETWEEN X(J) AND X(J+1) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) !DIMENSION X(1),Y(1),B(1),C(1),D(1) integer*4 N,N1,J,JR c real*8 X(N),Y(N),B(N),C(N),D(N) real*8 X,Y,B,C,D DIMENSION X(1),Y(1),B(1),C(1),D(1) real*8 S,R N1=N-1 S=0.0 DO 10 J=1,N1 D(J)=X(J+1)-X(J) R=(Y(J+1)-Y(J))/D(J) C(J)=R-S 10 S=R S=0.0 R=0.0 C(1)=0.0 C(N)=0.0 DO 20 J=2,N1 C(J)=C(J)+R*C(J-1) B(J)=(X(J-1)-X(J+1))*2.0-R*S S=D(J) 20 R=S/B(J) DO 30 JR=N1,2,-1 30 C(JR)=(D(JR)*C(JR+1)-C(JR))/B(JR) DO 40 J=1,N1 S=D(J) R=C(J+1)-C(J) D(J)=R/S C(J)=3.0*C(J) 40 B(J)=(Y(J+1)-Y(J))/S-(C(J)+R)*S RETURN END SUBROUTINE BSPOL(S,N,X,Y,B,C,D,T) implicit none C BINARY SEARCH, X ASCENDING OR DESCENDING C IMPLICIT DOUBLE PRECISION (A-H,O-Z) !DIMENSION X(1),Y(1),B(1),C(1),D(1) integer*4 N,IDIR,MLB,MUB,IDIR1,ML,MU,MAV c real*8 X(N),Y(N),B(N),C(N),D(N),T,S,Q real*8 X,Y,B,C,D,T,S,Q DIMENSION X(1),Y(1),B(1),C(1),D(1) IF (X(1).GT.X(N)) GO TO 10 IDIR=0 MLB=0 MUB=N GO TO 20 10 IDIR=1 MLB=N MUB=0 20 IDIR1=IDIR-1 IF (S.GE.X(MUB+IDIR)) GO TO 60 IF (S.LE.X(MLB-IDIR1)) GO TO 70 ML=MLB MU=MUB GO TO 40 30 IF (IABS(MU-ML).LE.1) GO TO 80 40 MAV=(ML+MU)/2 IF (S.LT.X(MAV)) GO TO 50 ML=MAV GO TO 30 50 MU=MAV GO TO 30 60 MU=MUB+IDIR+IDIR1 GO TO 90 70 MU=MLB-IDIR-IDIR1 GO TO 90 80 MU=MU+IDIR1 90 Q=S-X(MU) T=((D(MU)*Q+C(MU))*Q+B(MU))*Q+Y(MU) RETURN END SUBROUTINE LCOF (NMAX,X,F,B,C,D) implicit none C 26 JAN 93. SIMPLE LINEAR INTERPOLATION C IMPLICIT DOUBLE PRECISION (A-H,O-Z) !DIMENSION X(1), F(1), B(1), C(1), D(1) integer*4 NMAX,N c REAL*8 X(NMAX), F(NMAX), B(NMAX), C(NMAX), D(NMAX) REAL*8 X, F, B, C, D DIMENSION X(1), F(1), B(1), C(1), D(1) DO 10 N=1,NMAX-1 B(N)=(F(N+1)-F(N))/(X(N+1)-X(N)) C(N)=0.0 D(N)=0.0 10 CONTINUE RETURN END SUBROUTINE INTEG (DELTA,G,N,RESULT) implicit none C INCLUDES N=1 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer*4 N,NL1,NL2,K real*8 DELTA,result,SIGMA,SUM4,SUM2,SIG6 real*8 G(*) c real*8 G c dimension G(1) NL1=N-1 NL2=N-2 IF (REAL (N) -2.0*REAL (N/2)) 100,100,10 10 IF (N-1) 15,15,20 15 SIGMA=0.0 GO TO 70 20 IF(N-3) 30,30,40 30 SIGMA=G(1)+4.0*G(2)+G(3) GO TO 70 40 SUM4=0.0 DO 50 K=2,NL1,2 50 SUM4=SUM4+G(K) SUM2=0.0 DO 60 K=3,NL2,2 60 SUM2=SUM2+G(K) SIGMA=G(1)+4.0*SUM4+2.0*SUM2+G(N) 70 RESULT=DELTA*SIGMA RETURN 100 IF(N-2)110,110,120 110 SIGMA=1.5*(G(1)+G(2)) GO TO 70 120 IF(N-4)130,130,140 130 SIGMA=1.125*(G(1)+3.0*G(2)+3.0*G(3)+G(4)) GO TO 70 140 IF(N-6)150,150,160 150 SIGMA=G(1)+3.875*G(2)+2.625*G(3)+2.625*G(4)+3.875*G(5)+G(6) GO TO 70 160 IF (N-8)170,170,180 170 SIGMA=G(1)+3.875*G(2)+2.625*G(3)+2.625*G(4)+3.875*G(5)+2.0*G(6) 1 +4.0*G(7)+G(8) GO TO 70 180 SIG6=G(1)+3.875*G(2)+2.625*G(3)+2.625*G(4)+3.875*G(5)+G(6) SUM4=0.0 DO 190 K=7,NL1,2 190 SUM4=SUM4+G(K) SUM2=0.0 DO 200 K=8,NL2,2 200 SUM2=SUM2+G(K) SIGMA=SIG6+G(6)+4.0*SUM4+2.0*SUM2+G(N) GO TO 70 END