1 : C GAMTR SOURCE CHAT 05/01/13 00:15:48 5004 2 : SUBROUTINE GAMTR(S,DS,F1ST,F2ST,R1,R2,CC,SS,CS,ITR,IRZ,GAMMA, 3 : . PREC,RFSG,RFEP,RFPR,KERRE) 4 : C 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : -INC CCOPTIO 8 : C 9 : DIMENSION S(*),DS(*) 10 : DIMENSION CP(3) 11 : C 12 : IT=ITR+1 13 : GO TO(100,200),IT 14 : WRITE(IOIMP,20) IT 15 : KERRE=640 16 : RETURN 17 : C 18 : 100 CP(1)=S(1) 19 : CP(2)=S(2) 20 : CP(3)=S(4) 20 : CALL DIAGOD(CP) 21 : X=DS(4)*DS(4)-DS(1)*DS(2) 22 : R=MIN(R1,R2) 23 : Y=(R-S(1))*DS(2)+(R-S(2))*DS(1)+2.D0*S(4)*DS(4) 24 : Z=S(4)*S(4)-S(1)*S(2)+R*(S(1)+S(2))-R*R 25 : IF(ABS(X).LT.RFSG*RFSG) X=0.D0 26 : IF(X.EQ.0.D0) GO TO 1 27 : DELTA=Y*Y-4.D0*X*Z 28 : VRF=MAX(ABS(Y),RFSG) 29 : RFRF=VRF*VRF*RFPR*RFPR 30 : IF(IIMPI.EQ.9) WRITE(IOIMP,1001) X,Y,Z,DELTA 31 : IF(ABS(DELTA).LE.RFRF.AND.DELTA.LE.0.D0) DELTA=0.D0 32 : IF(DELTA.GE.0.D0) GO TO 2 33 : WRITE(IOIMP,22)DELTA 34 : GAMMA=100.D0 35 : RETURN 36 : 2 RADEL=DSQRT(DELTA) 37 : GAMT1=(-Y+RADEL)/(2.D0*X) 38 : GAMT2=(-Y-RADEL)/(2.D0*X) 39 : IF(IIMPI.EQ.9) WRITE(IOIMP,1002) GAMT1,GAMT2 40 : IF(GAMT1.GT.1.D0) GAMT1=1.D0 41 : IF(GAMT2.GT.1.D0) GAMT2=1.D0 42 : IF(GAMT1.EQ.1.D0.AND.GAMT2.LT.0.D0) GAMT2=0.D0 43 : IF(GAMT2.EQ.1.D0.AND.GAMT1.LT.0.D0) GAMT1=0.D0 44 : IF(GAMT1.GE.0.D0.OR.GAMT2.GE.0.D0) GO TO 8 45 : GAMT1=0.D0 46 : GAMT2=0.D0 47 : GO TO 8 48 : C 49 : 200 IF(F1ST.GT.R1) GO TO 3 50 : GAMT1=100.D0 51 : GO TO 4 52 : 3 F1DS=ROTA(DS,CC,SS,CS,1) 53 : IF(F1DS.NE.0.D0) GO TO 5 54 : GAMT1=100.D0 55 : GO TO 4 56 : 5 F1S=ROTA(S,CC,SS,CS,1) 57 : GAMT1=(R1-F1S)/F1DS 58 : IF(F1ST.GT.R1.AND.GAMT1.LT.0.D0) GAMT1=0.D0 59 : 4 IF(F2ST.GT.R2) GO TO 6 60 : GAMT2=100.D0 61 : GO TO 8 62 : 6 F2DS=ROTA(DS,CC,SS,CS,2) 63 : IF(F2DS.NE.0.D0) GO TO 7 64 : GAMT2=100.D0 65 : GO TO 8 66 : 7 F2S=ROTA(S,CC,SS,CS,2) 67 : GAMT2=(R2-F2S)/F2DS 68 : IF(F2ST.GT.R2.AND.GAMT2.LT.0.D0) GAMT2=0.D0 69 : C 70 : 8 GAGA=GAMT1*GAMT2 71 : IF(GAGA.LE.0.D0) GO TO 9 72 : GAMMA=MIN(GAMT1,GAMT2) 73 : IRZ=1 74 : IF(GAMMA.EQ.GAMT2) IRZ=2 75 : DENOR=MIN(GAMT1,GAMT2) 76 : DENOR=MAX(DENOR,RFPR) 77 : DIF=ABS(GAMT1-GAMT2)/DENOR 78 : IF(DIF.LE.PREC) IRZ=3 79 : RETURN 80 : 9 IF(GAMT1.EQ.0.D0.AND.GAMT2.EQ.0.D0) GO TO 12 81 : IF(GAMT1.LT.0.D0)GO TO 10 82 : GAMMA=GAMT1 83 : IRZ=1 84 : RETURN 85 : 10 GAMMA=GAMT2 86 : IRZ=2 87 : RETURN 88 : 12 GAMMA=0.D0 89 : IRZ=3 90 : RETURN 91 : 1 IF(Y.NE.0.D0) GO TO 11 92 : WRITE(IOIMP,21) 93 : GAMMA=100.D0 94 : RETURN 95 : 11 GAMMA=-Z/Y 96 : IRZ=1 97 : C 98 : 20 FORMAT(1X,'ERREUR DANS GAMTR IL FAUT DETERMINER IT =',I4) 99 : 21 FORMAT(1X,'ERREUR DANS GAMTR DETERMINANT EST NUL') 100 : 22 FORMAT(1X,'ERREUR DANS GAMTR DELTA EST NEGATIF DELTA =',1PE12.5) 101 : 1001 FORMAT(1X,'X =',1PD12.5,1X,'Y =',1PD12.5, 102 : . 1X,'Z =',1PD12.5,1X,'DELTA =',1PD12.5) 103 : 1002 FORMAT(1X,'GAMT1 =',1PD12.5,1X,'GAMT2 =',1PD12.5) 104 : C 105 : RETURN 106 : END 107 :
© Cast3M 2003 - All rights reserved.
Disclaimer