C--------1---------2---------3---------4---------5---------6---------7-- C C This program, PHOTONRAY.FOR, traces photon through the TIDE mirror C using Snell's law and the three best-fit parabolas of the mirror C surfaces. MIRRORFIT was used to calculate the parabolas. C C The program was written in Sept.1989 by Peggy Sloan. C C The program was modified in Jan.1990 to trace photons through TIDE C using the final Mirror configuration. C C The program was modified in Nov.90 to trace photons through the C flight version of TIDE. C C The program was modified in Dec.90 to use probability (random numbers) C to decide if and when a photon stops. For the Mirror and RPA grids, C the probability that the photon hits a wire is 10%, the probability C that a photon is reflected by any surface is 15%. C--------1---------2---------3---------4---------5---------6---------7-- DIMENSION Y0(3),X0(3),A0(3),X(50),Y(50),REF(4,3) CHARACTER*40 FNAME CHARACTER*1 NY CHARACTER*9 IDATE CHARACTER*8 ITIME C--------1---------2---------3---------4---------5---------6---------7-- C DATA REF/0.104,0.106,0.119,0.245,0.081,0.096,0.137,0.310, * 0.104,0.107,0.130,0.240/ C DATA REF/0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8/ RAD=180.0/3.14159268 YMN=5.03 YMX=7.4 PAST=0.10 PAST1=0.10 ITOF=0 IDEF=0 C RANDOM=SECNDS(0.0) IRAND=RANDOM C C SET UP CURVE COEFFICIENT ARRAYS (as calculated by MIRRORFIT) C A0(1)=1.14635563 X0(1)=3.00304532 Y0(1)=4.08808613 C A0(2)=1.19425845 X0(2)=2.82150245 Y0(2)=4.12670469 C A0(3)=1.26642263 X0(3)=2.57816219 Y0(3)=4.16414118 C C DETERMINE EQUATIONS ALONG MIRROR FRAMES C SLF=(5.03-4.869)/(3.3317-3.267) ! LOWER FRONT MIRROR FRAME BLF=5.03-SLF*3.3317 SUF=(7.5150-7.4)/(5.7930-5.6187) ! UPPER FONT MIRROR FRAME BUF=7.4-SUF*5.6187 C C OPEN THE PLOT C CALL OPNGKS 150 CALL DATE(IDATE) CALL TIME(ITIME) C C READ THE PHOTON STARTING RANGES C WRITE(6,210) 210 FORMAT(/'$Enter the name of the range file => ') READ(5,220) FNAME 220 FORMAT(A40) OPEN(UNIT=7,FILE=FNAME,TYPE='OLD',READONLY) READ(7,230) 230 FORMAT(1X) READ(7,240) ZI,ZMN,ZMX,AI,AMN,AMX 240 FORMAT(3G10.0) CLOSE(UNIT=7) IZ=ZI IA=AI C NPTS=IZ*IA C--------1---------2---------3---------4---------5---------6---------7-- C C SET UP PLOT C CALL SET (0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,1) CALL AGSETI ('FRAME.',2) CALL AGSETI ('SET.',-1) C CALL AGSETF ('GRID/TOP.',.89) CALL AGSETF ('GRID/BOTTOM.',.07) CALL AGSETF ('GRID/LEFT.',.07) CALL AGSETF ('GRID/RIGHT.',.89) C CALL AGSETF ('X/MIN.',0.00) CALL AGSETF ('X/MAX.',7.78) CALL AGSETF ('Y/MIN.',0.00) CALL AGSETF ('Y/MAX.',7.78) C CALL AGSETI('X/NICE.',0) CALL AGSETI('Y/NICE.',0) CALL ANOTAT('inches','inches',0,0,0,0) C--------1---------2---------3---------4---------5---------6---------7-- C X1=4.0 Y1=4.0 CALL EZXY (X1,Y1,1, * 'TIDE Flight Version - Photon Raytraces') C C READ AND PLOT TIDE MIRROR/RPA/TOF C CALL SETUSV('LW',1000) CALL PLREGB C--------1---------2---------3---------4---------5---------6---------7-- C C CALCULATE THE PHOTON RAY PATH C C--------1---------2---------3---------4---------5---------6---------7-- C SET UP INITIAL PHOTON PATH C C ICNT=0 ZINC=(ZMX-ZMN)/(IZ-1) AINC=(AMX-AMN)/(IA-1) DO 990 K=1,IZ DO 980 J=1,IA ICNT=ICNT+1 XST=6.0 YST=ZMN+(K-1)*ZINC ANGLE=AMN+(J-1)*AINC SST=TAN(-ANGLE/RAD) BST=YST-SST*XST C DO 930 I=1,NPTS X1=XST Y1=YST SNOW=SST BNOW=BST X(1)=X1 Y(1)=Y1 NR=1 C C DETERMINE IF THE PHOTON HITS THE LOWER MIRROR FRAME C 270 IF(NR.GT.3.AND.ABS(SNOW).GT.0.00001) THEN XTST=(4.66-BNOW)/SNOW IF(XTST.GT.4.543) THEN NR=NR+1 X(NR)=4.543 Y(NR)=4.543*SNOW+BNOW GO TO 920 ENDIF ENDIF C X2=(BLF-BNOW)/(SNOW-SLF) Y2=SLF*X2+BLF IFRAME=0 IF(Y2.LT.5.03.AND.Y2.GT.4.66.AND.X2.LT.3.3317) THEN IF(Y2.LT.4.8690) THEN X2=3.2670 Y2=SNOW*3.2670+BNOW THETA=2.0*ATAN(-SNOW) ELSE SM2=-1.0/SLF TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) ENDIF NR=NR+1 X(NR)=X2 Y(NR)=Y2 C T1=THETA/2.0*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ENDIF SNEW=(TAN(THETA)+SNOW)/(1.0-TAN(THETA)*SNOW) SNOW=SNEW BNOW=Y2-X2*SNOW IFRAME=1 ENDIF C C DETERMINE IF THE PHOTON HITS THE UPPER MIRROR FRAME C IF(IFRAME.EQ.0) THEN X2=(BUF-BNOW)/(SNOW-SUF) Y2=SUF*X2+BUF IF(Y2.GT.7.4.AND.X2.GT.4.6497.AND.X2.LT.6.055) THEN IF(Y2.GT.7.5150) THEN X2=5.7930 Y2=SNOW*5.7930+BNOW THETA=2.0*ATAN(SNOW) ELSE SM2=-1.0/SUF TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) ENDIF NR=NR+1 X(NR)=X2 Y(NR)=Y2 C T1=THETA/2.0*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ENDIF SNEW=(TAN(THETA)+SNOW)/(1.0-TAN(THETA)*SNOW) SNOW=SNEW BNOW=Y2-X2*SNOW IFRAME=1 ELSEIF(X2.GT.6.055)THEN NR=NR+1 X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW GO TO 920 ENDIF ENDIF IF(IFRAME.EQ.1) GO TO 400 C--------1---------2---------3---------4---------5---------6---------7-- C C CALCULATE THE INTERSECTION OF THE INITIAL PATH WITH THE MIRROR FRONT C IRMIR=1 C--------1---------2---------3---------4---------5---------6---------7-- 310 CALL MIRRORR(ANGLE,X(NR),Y(NR),X2,Y2,A0(1),X0(1),Y0(1), * SNOW,BNOW,YMN,YMX,SNEW,BNEW,IJUNC,THETA) IF(IJUNC.EQ.0) GO TO 920 NR=NR+1 IF(Y2.GT.7.4) THEN Y(NR)=7.4 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(7.4-BNOW)/SNOW ELSE NR=NR-1 ENDIF IF(X(NR).GT.6.055) THEN X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW ENDIF GO TO 920 ELSEIF(Y2.LT.5.03) THEN Y(NR)=5.03 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(5.03-BNOW)/SNOW ELSE NR=NR-1 ENDIF GO TO 920 ENDIF X(NR)=X2 Y(NR)=Y2 IF(RAN(IRAND).LT.PAST1) THEN SNOW=SNEW BNOW=BNEW C T1=THETA/2.0*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT.OR.NR.GT.30) THEN GO TO 920 ELSE IRMIR=-IRMIR IF(IRMIR.LT.0) GO TO 400 ENDIF ELSEIF(IRMIR.LT.0) THEN IRMIR=1 GO TO 400 ENDIF C C CALCULATE THE INTERSECTION OF THE INITIAL PATH WITH THE MIRROR MIDDLE C 320 CALL MIRRORR(ANGLE,X(NR),Y(NR),X2,Y2,A0(2),X0(2),Y0(2), * SNOW,BNOW,YMN,YMX,SNEW,BNEW,IJUNC,THETA) IF(IJUNC.EQ.0) GO TO 920 NR=NR+1 IF(Y2.GT.7.4) THEN Y(NR)=7.4 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(7.4-BNOW)/SNOW ELSE NR=NR-1 ENDIF IF(X(NR).GT.6.055) THEN X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW ENDIF GO TO 920 ELSEIF(Y2.LT.5.03) THEN Y(NR)=5.03 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(5.03-BNOW)/SNOW ELSE NR=NR-1 ENDIF GO TO 920 ENDIF X(NR)=X2 Y(NR)=Y2 IF(RAN(IRAND).LT.PAST1) THEN SNOW=SNEW BNOW=BNEW C T1=THETA/2.0*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT.OR.NR.GT.30) THEN GO TO 920 ENDIF IRMIR=-IRMIR IF(IRMIR.LT.0) GO TO 310 ELSEIF(IRMIR.LT.0) THEN GO TO 310 ENDIF C C CALCULATE THE INTERSECTION OF THE INITIAL PATH WITH THE MIRROR BACK C 330 CALL MIRRORR(ANGLE,X(NR),Y(NR),X2,Y2,A0(3),X0(3),Y0(3), * SNOW,BNOW,YMN,YMX,SNEW,BNEW,IJUNC,THETA) IF(IJUNC.EQ.0) GO TO 920 NR=NR+1 IF(Y2.GT.7.4) THEN Y(NR)=7.4 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(7.4-BNOW)/SNOW ELSE NR=NR-1 ENDIF IF(X(NR).GT.6.055) THEN X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW ENDIF GO TO 920 ELSEIF(Y2.LT.5.03) THEN Y(NR)=5.03 IF(ABS(SNOW).GT.0.00001) THEN X(NR)=(5.03-BNOW)/SNOW ELSE NR=NR-1 ENDIF GO TO 920 ENDIF X(NR)=X2 Y(NR)=Y2 IF(RAN(IRAND).LT.PAST1) THEN SNOW=SNEW BNOW=BNEW C T1=THETA/2.0*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT.OR.NR.GT.30) THEN GO TO 920 ENDIF IRMIR=-1 GO TO 320 ELSE NR=NR+1 X(NR)=2.33 Y(NR)=2.33*SNOW+BNOW IF(Y(NR).GT.7.78.OR.Y(NR).LT.4.0) THEN Y(NR)=7.78 IF(ABS(SNOW).LT.0.00001) THEN NR=NR-1 ELSE X(NR)=(7.78-BNOW)/SNOW ENDIF ENDIF GO TO 920 ENDIF IF(NR.GT.30) THEN GO TO 920 ENDIF C C DID PHOTON ENTER RPA? C 400 IF(Y(NR).LT.5.5.AND.Y(NR).GT.5.03) THEN YTST=3.3317*SNOW+BNOW IF(YTST.LT.5.03) THEN ! PHOTON HITS LOWER FRAME NR=NR+1 Y(NR)=5.03 IF(ABS(SNOW).LT.0.00001) THEN NR=NR-1 ELSE X(NR)=(5.03-BNOW)/SNOW ENDIF GO TO 920 ENDIF ENDIF NR=NR+1 Y(NR)=4.66 IF(ABS(SNOW).LT.0.00001) THEN NR=NR-1 ELSE X(NR)=(4.66-BNOW)/SNOW ENDIF IRMIR=1 IRPA=1 C IF(X(NR).GT.6.055.OR.X(NR).LT.2.33) THEN ! PHOTON EXITS AT APERTURE X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW GO TO 920 ELSEIF(X(NR).LT.3.317) THEN X(NR)=3.317 Y(NR)=3.317*SNOW+BNOW GO TO 920 ELSEIF(X(NR).GT.4.543) THEN ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT.OR.NR.GT.30) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) THEN XTST=5.6187 ELSE XTST=(7.4-BNOW)/SNOW ENDIF IF(XTST.LT.5.6187) THEN GO TO 310 ELSE NR=NR+1 X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW GO TO 920 ENDIF ENDIF ENDIF YTST=4.543*SNOW+BNOW IF(YTST.LT.4.66.AND.YTST.GT.4.56) THEN NR=NR+1 Y(NR)=YTST X(NR)=4.543 GO TO 920 ENDIF C--------1---------2---------3---------4---------5---------6---------7-- C C TRACE PHOTON THROUGH THE RPA C C--------1---------2---------3---------4---------5---------6---------7-- C C FIRST GRID C IRPA=1 410 IF(RAN(IRAND).GT.PAST) THEN IF(IRPA.GT.0) THEN GO TO 420 ELSE GO TO 270 ENDIF ELSE IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=4.56 XTST=(4.56-BNOW)/SNOW X(NR)=XTST C--------1---------2---------3---------4---------5---------6---------7-- ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF((XTST.GT.4.543.OR.XTST.LT.3.317).OR.(RAN1.GT.REFLECT) * .OR.(NR.GT.30)) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(IRPA.GT.0) THEN IF(ABS(SNOW).LT.0.00001) THEN XTST=5.1687 ELSE XTST=(7.4-BNOW)/SNOW ENDIF IF(XTST.LT.5.6187) THEN IRMIR=1 GO TO 270 ELSE NR=NR+1 X(NR)=6.055 Y(NR)=6.055*SNOW+BNOW GO TO 920 ENDIF ELSE IRPA=-IRPA ENDIF ENDIF ENDIF C C SECOND GRID C 420 IF(RAN(IRAND).GT.PAST) THEN IF(IRPA.GT.0) THEN GO TO 430 ELSE GO TO 410 ENDIF ELSE IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=4.5288 XTST=(4.5288-BNOW)/SNOW X(NR)=XTST ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF((XTST.GT.4.543.OR.XTST.LT.3.317).OR.(RAN1.GT.REFLECT) * .OR.(NR.GT.30)) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IRPA=-IRPA IF(IRPA.LT.0) GO TO 410 ENDIF ENDIF C C THIRD GRID C 430 IF(RAN(IRAND).GT.PAST) THEN IF(IRPA.GT.0) THEN GO TO 440 ELSE GO TO 270 ENDIF ELSE IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=4.3418 XTST=(4.3418-BNOW)/SNOW X(NR)=XTST ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF((XTST.GT.4.543.OR.XTST.LT.3.317).OR.(RAN1.GT.REFLECT) * .OR.(NR.GT.30)) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IRPA=-IRPA IF(IRPA.LT.0) GO TO 420 ENDIF ENDIF C C FOURTH GRID C 440 IF(RAN(IRAND).GT.PAST) THEN IF(IRPA.GT.0) THEN GO TO 460 ELSE GO TO 270 ENDIF ELSE IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=4.3105 XTST=(4.3105-BNOW)/SNOW X(NR)=XTST ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF((XTST.GT.4.543.OR.XTST.LT.3.317).OR.(RAN1.GT.REFLECT) * .OR.(NR.GT.30)) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IRPA=-IRPA IF(IRPA.LT.0) GO TO 430 ENDIF ENDIF C C DID PHOTON ENTER DEFLECTOR PORTION OF TOF? C 460 IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=4.0 X(NR)=(4.0-BNOW)/SNOW IF(X(NR).LT.3.317.OR.X(NR).GT.4.543) GO TO 920 ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 IF(X(NR).LT.3.733.OR.X(NR).GT.4.127) THEN REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IRPA=-1 IF(ABS(SNOW).LT.0.00001) GO TO 920 XTST=(4.3105-BNOW)/SNOW IF(XTST.LT.3.317) THEN NR=NR+1 X(NR)=3.317 Y(NR)=3.317*SNOW+BNOW GO TO 920 ELSEIF(XTST.GT.4.543)THEN NR=NR+1 X(NR)=4.543 Y(NR)=4.543*SNOW+BNOW GO TO 920 ELSE GO TO 440 ENDIF ENDIF ELSE IF(RAN(IRAND).LT.PAST) THEN REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IRPA=-1 GO TO 440 ENDIF ENDIF ENDIF C C TRACE PHOTON THROUGH THE LOWER IMMERSION LENS C IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=3.7244 X(NR)=(3.7244-BNOW)/SNOW IF(X(NR).LT.3.733.OR.X(NR).GT.4.127) THEN ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 X(NR)=(3.9375-BNOW)/SNOW Y(NR)=3.9375 GO TO 920 ENDIF ENDIF C--------1---------2---------3---------4---------5---------6---------7-- C C TRACE PHOTON THROUGH THE CURVED DEFLECTORS C C--------1---------2---------3---------4---------5---------6---------7-- C DID PHOTON GET PASSED THE TOP INNER DEFLECTION LENS? C NR=NR+1 IF(SNOW.GT.0.0) THEN Y(NR)=3.4488 X(NR)=(3.4488-BNOW)/SNOW IF(X(NR).GT.3.6560) THEN GO TO 520 ELSE ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 X(NR)=(3.6619-BNOW)/SNOW Y(NR)=3.6619 GO TO 920 ENDIF ENDIF ELSE IF(ABS(SNOW).LT.0.00001) GO TO 920 Y(NR)=3.4173 X(NR)=(3.4173-BNOW)/SNOW IF(X(NR).GT.4.2008) THEN ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 X(NR)=(3.6619-BNOW)/SNOW Y(NR)=3.6619 GO TO 920 ENDIF ENDIF ENDIF 500 CALL ARCINT(3.5270,3.450,0.129,SNOW,BNOW,X1,Y1,X2,Y2,ISOLV) IF(ISOLV.EQ.1) THEN IF(X1.GT.3.5270.AND.X1.LT.3.656.AND. * Y1.LT.3.4488.AND.Y1.GT.3.321) THEN NR=NR+1 X(NR)=X1 Y(NR)=Y1 GO TO 920 ELSEIF(X2.GT.3.5270.AND.X2.LT.3.656.AND. * Y2.LT.3.4488.AND.Y2.GT.3.321) THEN NR=NR+1 X(NR)=X2 Y(NR)=Y2 GO TO 920 ENDIF ENDIF C C DID PHOTON GET PASSED THE TOP OUTER DEFLECTION LENS? C 520 Y4=3.4173 IF(ABS(SNOW).LT.0.00001) GO TO 920 X4=(Y4-BNOW)/SNOW IF(X4.GT.4.2008) THEN NR=NR+1 X(NR)=X4 Y(NR)=Y4 GO TO 920 ENDIF 530 CALL ARCINT(3.5650,3.4173,0.63583,SNOW,BNOW,X1,Y1,X2,Y2,ISOLV) IHIT=0 IF(ISOLV.EQ.1) THEN IF(X1.GT.3.7484.AND.X1.LT.4.2008.AND. * Y1.LT.3.4173.AND.Y1.GT.2.7638) THEN NR=NR+1 X(NR)=X1 Y(NR)=Y1 IHIT=1 ELSEIF(X2.GT.3.7484.AND.X2.LT.4.2008.AND. * Y2.LT.3.4173.AND.Y2.GT.2.7638) THEN NR=NR+1 X(NR)=X2 Y(NR)=Y2 IHIT=1 ENDIF IF(IHIT.EQ.1) THEN SM2=(Y(NR)-3.4173)/(X(NR)-3.5650) ANG=ATAN((SM2-SNOW)/(1.0+SM2*SNOW)) T1=ANG*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE SM2=(3.4173-Y(NR))/(3.5650-X(NR)) TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) TANG1=TAN(THETA) SNEW=(TANG1+SNOW)/(1.0-TANG1*SNOW) SNOW=SNEW BNOW=Y(NR)-SNOW*X(NR) YTST=3.2854*SNOW+BNOW IF(YTST.GT.3.0079) GO TO 500 YTST=2.8575*SNOW+BNOW IF(YTST.GT.2.4508) GO TO 550 IF(ABS(SNOW).LT.0.00001) GO TO 920 XTST=(2.7638-BNOW)/SNOW IF(XTST.GT.3.7484) GO TO 530 ENDIF ENDIF ENDIF C C DID PHOTON GET PASSED THE BOTTOM OUTER DEFLECTION LENS? C Y4=2.4818 IF(ABS(SNOW).LT.0.00001) GO TO 920 X4=(Y4-BNOW)/SNOW IF(X4.GT.3.5848) THEN NR=NR+1 X(NR)=X4 Y(NR)=Y4 GO TO 920 ENDIF CALL ARCINT(3.472,2.418,0.129,SNOW,BNOW,X1,Y1,X2,Y2,ISOLV) IHIT=0 IF(ISOLV.EQ.1) THEN IF(X1.GT.3.3437.AND.X1.LT.3.5848.AND. * Y1.LT.2.547.AND.Y1.GT.2.4193) THEN NR=NR+1 X(NR)=X1 Y(NR)=Y1 IHIT=1 ELSEIF(X2.GT.3.3437.AND.X2.LT.3.5848.AND. * Y2.LT.2.547.AND.Y2.GT.2.4193) THEN NR=NR+1 X(NR)=X2 Y(NR)=Y2 IHIT=1 ENDIF IF(IHIT.EQ.1) THEN SM2=(Y(NR)-2.418)/(X(NR)-3.4720) ANG=ATAN((SM2-SNOW)/(1.0+SM2*SNOW)) T1=ANG*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE SM2=(2.418-Y(NR))/(3.472-X(NR)) TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) TANG1=TAN(THETA) SNEW=(TANG1+SNOW)/(1.0-TANG1*SNOW) SNOW=SNEW BNOW=Y(NR)-SNOW*X(NR) IF(ABS(SNOW).LT.0.00001) GO TO 920 X1=(3.0079-BNOW)/SNOW IF(X1.GT.3.656) THEN NR=NR+1 X(NR)=(2.7638-BNOW)/SNOW Y(NR)=2.7638 GO TO 920 ELSEIF(X1.GT.3.2854) THEN NR=NR+1 X(NR)=(3.3863-BNOW)/SNOW Y(NR)=3.3863 GO TO 920 ELSE GO TO 550 ENDIF ENDIF ENDIF ENDIF C C DID PHOTON GET PASSED THE BOTTOM INNER DEFLECTION LENS? C Y4=3.0472 IF(ABS(SNOW).LT.0.00001) GO TO 920 X4=(Y4-BNOW)/SNOW IF(X4.LT.3.2665) THEN NR=NR+1 X(NR)=X4 Y(NR)=Y4 GO TO 920 ENDIF C S2=(3.0472-3.0079)/(3.2665-3.2854) B2=(3.0472-3.2665) X4=(B2-BNOW)/(SNOW-S2) IF(X4.LT.3.2854.AND.X4.GT.3.2665) THEN NR=NR+1 X(NR)=X4 Y(NR)=SNOW*X4+BNOW GO TO 920 ENDIF C 550 CALL ARCINT(3.4343,2.4508,0.57677,SNOW,BNOW,X1,Y1,X2,Y2,ISOLV) IHIT=0 IF(ISOLV.EQ.1) THEN IF(X1.GT.2.8575.AND.X1.LT.3.2854.AND. * Y1.LT.3.0079.AND.Y1.GT.2.4508) THEN NR=NR+1 X(NR)=X1 Y(NR)=Y1 IHIT=1 ELSEIF(X2.GT.2.8575.AND.X2.LT.3.2854.AND. * Y2.LT.3.0079.AND.Y2.GT.2.4508) THEN NR=NR+1 X(NR)=X2 Y(NR)=Y2 IHIT=1 ENDIF IF(IHIT.EQ.1) THEN SM2=(Y(NR)-2.4508)/(X(NR)-3.4343) ANG=ATAN((SM2-SNOW)/(1.0+SM2*SNOW)) T1=ANG*RAD ANG=ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE SM2=(2.4508-Y(NR))/(3.4343-X(NR)) TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) TANG1=TAN(THETA) SNEW=(TANG1+SNOW)/(1.0-TANG1*SNOW) SNOW=SNEW BNOW=Y(NR)-SNOW*X(NR) IF(ABS(SNOW).LT.0.00001) GO TO 920 X1=(2.1575-BNOW)/SNOW IF(X1.LT.2.925.OR.X1.GT.3.325) GO TO 920 ENDIF ENDIF ENDIF C--------1---------2---------3---------4---------5---------6---------7-- C C DID PHOTON ENTER TOF? C IF(ABS(SNOW).LT.0.00001) GO TO 920 NR=NR+1 Y(NR)=2.1575 X(NR)=(2.1575-BNOW)/SNOW IF(X(NR).LT.2.925.OR.X(NR).GT.3.325) THEN ANG=ATAN(SNOW) T1=ANG*RAD ANG=90.0-ABS(T1)-0.5 IPER=MOD(ANG,15.0)+1 IF(IPER.GT.4) IPER=4 REFLECT=REF(IPER,2) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) GO TO 920 IF(SNOW.LT.0.0) THEN NR=NR+1 Y(NR)=2.4193 X(NR)=(2.4193-BNOW)/SNOW GO TO 920 ELSE NR=NR+1 Y(NR)=2.4508 X(NR)=(2.4508-BNOW)/SNOW GO TO 920 ENDIF ENDIF ELSE IF(RAN(IRAND).LT.PAST) THEN REFLECT=REF(IPER,1) RAN1=RAN(IRAND) IF(RAN1.GT.REFLECT) THEN GO TO 920 ELSE XDIF=X(NR)-X(NR-1) X2=X(NR)+XDIF Y2=Y(NR-1) SNOW=(Y2-Y(NR))/XDIF BNOW=Y2-SNOW*X2 IF(ABS(SNOW).LT.0.00001) GO TO 920 IF(SNOW.LT.0.0) THEN NR=NR+1 Y(NR)=2.4193 X(NR)=(2.4193-BNOW)/SNOW GO TO 920 ELSE NR=NR+1 Y(NR)=2.4508 X(NR)=(2.4508-BNOW)/SNOW GO TO 920 ENDIF ENDIF ENDIF ENDIF C C TRACE PHOTON THROUGH THE TOF C Y4=0.150 IF(ABS(SNOW).LT.0.00001) GO TO 920 X4=(Y4-BNOW)/SNOW IF(X4.GT.3.449) THEN Y5=0.8091 X5=(Y5-BNOW)/SNOW IF(X5.GT.3.449) THEN X(NR)=X5 Y(NR)=Y5 GO TO 920 ENDIF ENDIF C IF(X4.LT.2.488) THEN X4=2.488 Y4=SNOW*X4+BNOW X(NR)=X4 Y(NR)=Y4 GO TO 920 ENDIF C X(NR)=X4 Y(NR)=Y4 C C--------1---------2---------3---------4---------5---------6---------7-- 920 IF(Y(NR).LT.4.0) THEN IDEF=IDEF+1 IF(Y(NR).LT.2.1575) THEN CALL CURVE(X,Y,NR) WRITE(55,955) ICNT,X(1),Y(1),ANGLE,SST,BST,NR,X(NR),Y(NR) 955 FORMAT(' I,X,Y,A,M,B=',I8,5F12.6/' N,X,Y=',I8,2F12.6) ITOF=ITOF+1 ENDIF ENDIF 980 CONTINUE 990 CONTINUE WRITE(55,956) ICNT,IDEF,ITOF 956 FORMAT(I8,' photons started'/I8,' photons made it through the ', *'mirror and RPA'/I8,' photons entered the TOF') C C--------1---------2---------3---------4---------5---------6---------7-- C FINISH PLOT C CALL SET (0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,1) CALL WTSTR(0.90,0.915,IDATE,6,0,0) CALL WTSTR(0.90,0.900,ITIME,6,0,0) CALL FRAME C C CREATE NEXT PLOT OR STOP EXECUTION C WRITE(6,800) 800 FORMAT(/'$Do you want another plot run? (Y/N) => ') READ(5,810) NY 810 FORMAT(A1) IF(NY.EQ.'Y'.OR.NY.EQ.'y') GO TO 150 C CALL CLSGKS CLOSE(UNIT=8) C STOP END C--------*---------*---------*---------*---------*---------*---------*-- C--------*---------*---------*---------*---------*---------*---------*-- C C SUBROUTINE: ARCINT C C This subroutine calculates the 2 points where a line intersects a C circle, given the center and radius of the circle and the slope and C Y-intercept of the line. C--------*---------*---------*---------*---------*---------*---------*-- SUBROUTINE ARCINT(CX,CY,RADIUS,SLOPE,YINT,X1,Y1,X2,Y2,ISOLV) C C SET UP QUADRATIC EQUATIONS COEFFICIENTS FOR THE FOLLOWING EQUATION C C (X-CX)**2 + (Y-CY)**2 = RADIUS WITH Y = SLOPE*X +YINT C A=SLOPE*SLOPE+1.0 BCY=YINT-CY B=(SLOPE*BCY-CX)*2.0 C=BCY*BCY+CX*CX-RADIUS*RADIUS C C SOLVE QUADRADIC EQUATION FOR BOTH X INTERCEPT POINTS C A2=2.0*A DIST=(B*B-4.0*A*C) IF(DIST.LT.0.000001) THEN ISOLV=0 ELSE DIST=SQRT(B*B-4.0*A*C) X1=(-B-DIST)/A2 X2=(-B+DIST)/A2 Y1=SLOPE*X1+YINT Y2=SLOPE*X2+YINT ISOLV=1 ENDIF C RETURN END C--------*---------*---------*---------*---------*---------*---------*-- C--------*---------*---------*---------*---------*---------*---------*-- C C SUBROUTINE: MIRRORR C C This subroutine calculates the point and the line of a photon and a C TIDE mirror surface. C--------*---------*---------*---------*---------*---------*---------*-- SUBROUTINE MIRRORR(ANG,X1,Y1,X2,Y2,AM,XM,YM,SNOW,BNOW,YMN,YMX, * SNEW,BNEW,IJUNC,THETA) C YMN=4.66 YMX=7.77 IJUNC=1 IF(ABS(SNOW).LT.0.00001) THEN Y2=Y1 X2=(Y2-YM)**2/(4.0*AM)+XM ELSE B=-1.0*(2.0*YM+4.0*AM/SNOW) C=YM*YM-4.0*AM*(X1-Y1/SNOW-XM) ROOT=B*B-4.0*C IF(ROOT.LT.0.00001) THEN IJUNC=0 RETURN ELSE YP=(-B+SQRT(B*B-4.0*C))/2.0 YN=(-B-SQRT(B*B-4.0*C))/2.0 ITST=0 ENDIF C IF(YN.GT.YMN.AND.YN.LT.YMX) THEN Y2=YN ITST=ITST+1 ENDIF IF(YP.GT.YMN.AND.YP.LT.YMX) THEN Y2=YP ITST=ITST+1 ENDIF C IF(ITST.EQ.0) THEN ! NO INTERCEPT IJUNC=0 RETURN ELSEIF(ITST.GT.1) THEN D1=ABS(Y1-YP) D2=ABS(Y1-YN) IF(D1.LT.D2) THEN Y2=YP ELSE Y2=YN ENDIF ENDIF X2=(Y2-Y1)/SNOW+X1 ENDIF C C CALCULATE THE PERPENDICULAR TO THE TANGENT AT (X2,Y2) C IF(X2.LE.XM) THEN IJUNC=0 RETURN ELSE TANG=SQRT(AM/(X2-XM)) SM2=-1.0/TANG ENDIF C C CALCULATE THE SLOPE OF THE REFLECTED RAY PATH (SNELL'S LAW) C TANG=(SM2-SNOW)/(1.0+SM2*SNOW) THETA=2.0*ATAN(TANG) SNEW=(TAN(THETA)+SNOW)/(1.0-TAN(THETA)*SNOW) BNEW=Y2-X2*SNEW C--------1---------2---------3---------4---------5---------6---------7-- RETURN END C C--------*---------*---------*---------*---------*---------*---------*-- C--------*---------*---------*---------*---------*---------*---------*-- C C SUBROUTINE: PLREGB C C This subroutine is a combination of TIDEPLOT subroutines RDREGB and C PLREGB. They were modified by Peggy Sloan in Sept.1989 for use in C PHOTONRAY. C--------*---------*---------*---------*---------*---------*---------*-- SUBROUTINE PLREGB C CHARACTER*40 ANAME,ATIT CHARACTER*1 YN C DIMENSION XRB(50,50),YRB(50,50),NRB(50),XR(50),YR(50) C NAMELIST /REG/ NREG,NCELL,IREG,MAT,CUR,DEN,ITRI, X IBOUND,DX,XMIN,XMAX,DY,YMIN,YMAX,YREG1, X YREG2,NPOINT,IPRINT,XREG1,XREG2,KREG1,KREG2,KMAX, X LREG1,LREG2,LMAX,CONV,LINX,LINY,NDRIVE C NAMELIST /PO/ NT,X0,Y0,X,Y,R,THETA,NEW C C INPUT THE 'AUTOMESH' INPUT FILE NAME AND OPEN THE FILE C WRITE(6,110) 110 FORMAT(/' Enter the name of the AUTOMESH input file.'/ * '$(default: TIDE.DAT) ==> ') READ(5,120) LN,ANAME 120 FORMAT(Q,A40) IF(LN.LT.1) ANAME='TIDE.DAT' OPEN(UNIT=1,FILE=ANAME,TYPE='OLD') C C READ TITLE C READ(1,150) LN,ATIT 150 FORMAT(1X,Q,A40) C--------*---------*---------*---------*---------*---------*---------*-- C C READ THE REGIONAL BOUNDARY DATA POINTS C PI=3.1415927 RAD=PI/180.0 CONV=1.0 CUR=0.0 READ(1,REG) C C READ BOUNDARY CONDITIONS (NAMELIST REG) C IF(NREG.GT.50) NREG=50 CUR=0.0 DO I=1,NREG CUR=0.0 IF(I.GT.1) THEN READ(1,REG) ENDIF C C READ EACH POINT IN THE BOUNDARY C NPT=0 DO J=1,NPOINT C C SET $PO NAMELIST DEFAULTS C NT=1 X0=0.0 Y0=0.0 X=-99999.0 Y=-99999.0 R=-99999.0 THETA=-99999.0 NEW=0 C READ(1,PO) C C CALCULATE POINTS ALONG A STRAIGHT LINE SEGMENT OF THE BOUNDARY C IF(NT.EQ.1) THEN NPT=NPT+1 IF(R.LT.-90000.0) THEN XRB(I,NPT)=X+X0 YRB(I,NPT)=Y+Y0 ELSE TH=THETA*RAD XRB(I,NPT)=R*COS(TH)+X0 YRB(I,NPT)=R*SIN(TH)+Y0 ENDIF C C CALCULATE POINTS ALONG AN ARC SEGMENT OF THE BOUNDARY C ELSEIF(NT.EQ.2) THEN IF(NPT.EQ.0) THEN WRITE(6,225) I,J STOP ELSEIF(R.LT.-90000.0.OR.THETA.LT.-90000.0) THEN WRITE(6,210) I,J 210 FORMAT(/' The arc in region ',I2,', at point ',I2,' is ' * 'missing the R or THETA value.') STOP ELSE DELX=XRB(I,NPT)-X0 DELY=YRB(I,NPT)-Y0 IF(DELX.EQ.0.0.AND.DELY.EQ.0.0) THEN WRITE(6,220) I,J 220 FORMAT(/' There is a error in the location of (X0,Y0)', * ' for the arc'/' in region ',I2,' at point ',I2,'.') STOP ELSEIF(DELX.EQ.0.0) THEN IF(DELY.LT.0.0) THEN !CALCULATE THE R,THETA OF THE THST=-90.0*RAD !FIRST POINT ON THE ARC ELSE THST=90.0*RAD ENDIF ELSE THST=ATAN(DELY/DELX) ENDIF C FIVE=5.0*RAD DEG=THETA*RAD-THST ! CALCULATE THE ANGLE FORMED BY ARC JPT=ABS(DEG)/FIVE ! CALCULATE POINTS EVERY 5 DEG ON ARC IF(JPT*FIVE.LT.ABS(DEG)) JPT=JPT+1 DELTA=DEG/JPT DO L=1,JPT THST=THST+DELTA NPT=NPT+1 XRB(I,NPT)=R*COS(THST)+X0 YRB(I,NPT)=R*SIN(THST)+Y0 ENDDO ENDIF C C CALCULATE POINTS ALONG A HYPERBOLIC SEGMENT OF THE BOUNDARY C ELSEIF(NT.EQ.3) THEN IF(NPT.EQ.0) THEN WRITE(6,225) I,J 225 FORMAT(/' The first point in a region can not be an arc', * ' or a hyperbola.'/' Check the first point of region',I2) STOP ELSEIF(X0.NE.0.0.OR.Y0.NE.0.0.OR.THETA.GT.-90000.0) THEN WRITE(6,230) I,J 230 FORMAT(/' The only acceptable input for a hyperbola is ', * 'X, Y, R, and NT=3.'/' Check your hyperbola in region ', * I2,' at point ',I2,'.') STOP ELSEIF(R.LE.0.0.OR.X.EQ.0.0.OR.Y.EQ.0.0) THEN WRITE(6,240) I,J 240 FORMAT(/' You must specify the minimum distance from the', * ' origin to your hyperbola.'/' Neither X nor Y can be ', * 'xero.'/' Check region ',I2,' at point ',I2,'.') STOP ELSE YTST=R*R/2.0/X IF(ABS(YTST-Y).GT.0.01) THEN WRITE(6,250) I,J 250 FORMAT(/' Hyperbolas must obey the equation 2*X*Y=R*R'/ * ' Check region ',I2,' at point ',I2,'.') STOP ELSE X1=XRB(I,NPT) DELX=(X-X1)/20.0 XNUM=R*R/2.0 DO K=1,20 IF(K.EQ.20) THEN X1=X ELSE X1=X1+DELX ENDIF IF(X1.EQ.0.0) THEN WRITE(6,260) I,J 260 FORMAT(/' Your hyperbola can not cross an axis.'/ * ' Check region ',I2,' at point ',I2,'.') STOP ELSE NPT=NPT+1 YRB(I,NPT)=XNUM/X1 XRB(I,NPT)=X1 ENDIF ENDDO ENDIF ENDIF ENDIF C ENDDO NRB(I)=NPT ENDDO CLOSE(UNIT=1) C C PLOT THE REGIONAL BOUNDARIES C CALL DASHDB('73567'O) DO I=2,NREG NR=NRB(I) DO J=1,NR XR(J)=XRB(I,J) YR(J)=YRB(I,J) ENDDO CALL CURVED(XR,YR,NR) ENDDO C RETURN END