1 : C SIGELP SOURCE CHAT 05/01/13 03:16:35 5004 2 : SUBROUTINE SIGELP(SIGMA,DSIGMA,SIGEL,DSIGP,STOT,STEST,ITYPE, 3 : 1 SEL,SPLA,IBOU,S0,SELAS,XMAT,COVNMS,ALFAH,IMAPLA,SSTAR) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : DIMENSION SIGMA(*),DSIGMA(*),STOT(*),SIGEL(*),DSIGP(*) 7 : 1,XMAT(*),COVNMS(*) 8 : DATA UN/1.0D0/,ZERO/0.0D0/ 9 : C 10 : CCCCCCCCCCCCCCC CALCUL DE L INTERSECTION AVEC LE CONVEXE 11 : C 11 : CALL ZDANUL(SIGEL,6) 11 : CALL ZDANUL(DSIGP,6) 12 : AA=VONMIS(DSIGMA,ITYPE,ALFAH,COVNMS) 13 : BB=VNMS12(DSIGMA,SIGMA,ITYPE,ALFAH,COVNMS) 14 : AA=AA*AA 15 : AA1=AA 16 : BB1=BB 17 : TST=STEST*STEST 18 : IF(IMAPLA.EQ.11) GO TO 20 19 : C=S0*S0 20 : CC=C-SELAS*SELAS 21 : IF(IMAPLA.NE.5) GO TO 11 22 : C 23 : CCCCCCCCCCCCCC CAS DRUCKER PRAGER 24 : C 25 : EE=TRACE(DSIGMA) 26 : FF=TRACE(SIGMA) 27 : XMAT22=XMAT(2)*XMAT(2) 28 : XMAT11=XMAT(1)*XMAT(1) 29 : XMAT13=XMAT(1)*SELAS 30 : AA=AA*XMAT22-EE*EE*XMAT11 31 : BB=BB*XMAT22-EE*FF*XMAT11+XMAT13*EE 32 : GO TO 11 33 : 20 CONTINUE 34 : C 35 : CCCCCCCCCCCCAS CHEN ET CHEN 36 : C 37 : IN=0 38 : 90 CONTINUE 39 : IF (IN.EQ.0) GO TO 120 40 : AAA=AA 41 : BBB=BB 42 : CCC=CC 43 : DDD=DD 44 : BB=BB1 45 : AA=AA1 46 : 120 CONTINUE 47 : DINV1=TRACE(DSIGMA) 48 : AA=AA/3+IN*DINV1*DINV1/6 49 : CC=VONMIS(SIGMA,ITYPE,ALFAH,COVNMS) 50 : DINVO=TRACE(SIGMA) 51 : CC=CC*CC/3+IN*DINVO*DINVO/6+XMAT(2)*DINVO/3 52 : CC=CC+(-UN+XMAT(1)*DINVO/3)*SELAS*SELAS 53 : BB=BB*2/3+IN*DINV1*DINVO/3+DINV1*XMAT(2)/3 54 : BB=BB+XMAT(1)*SELAS*SELAS*DINV1/3 55 : BB=BB/2 56 : CCCCC 57 : 11 CONTINUE 58 : IF(ABS(AA).GT.TST) GO TO 1 59 : C 60 : CCCCCCCCCC AA EST NUL UNE SEULE SOLUTION 61 : C 62 : IF(ABS(BB).GT.TST) GO TO 10 63 : C 64 : CCCCCCCCCC BB EST NUL <=> DSIGMA=0 SIGMA0=STOT 65 : C 66 : 110 CONTINUE 67 : XX=SELAS/(S0+1.E-20) 67 : CALL AEQBX(SIGEL,SIGMA,XX,IBOU) 67 : CALL AEQBPC(DSIGP,SIGMA,SIGEL,UN,-UN,IBOU) 68 : SEL=0. 69 : SPLA=1. 70 : GO TO 555 71 : 10 CONTINUE 72 : C 73 : CCCCCCCCCCCCCCAS BB NON NUL 74 : C 75 : SEL=CC/(2.*BB) 76 : GO TO 333 76 : 333 CALL AEQBPC(SIGEL,SIGMA,DSIGMA,UN,SEL,IBOU) 77 : SPLA=1.-SEL 77 : CALL AEQBX(DSIGP,DSIGMA,SPLA,IBOU) 78 : GO TO 555 79 : 1 CONTINUE 80 : C 81 : CCCCCCCCCCCCCC AA NON NUL ON A 2 SOLUTIONS 82 : C 83 : IF(ABS(CC).GT.TST) GO TO 2 84 : C 85 : CCCCCCCCCCCCCC 1 SEULE SOLUTION CC NUL 86 : C 87 : SEL=0. 88 : GO TO 333 89 : 2 CONTINUE 90 : C 91 : CCCCCCCCCCCC 2 SOLUTIONS 92 : C 93 : DD=4.*BB*BB-4.*AA*CC 94 : B=2.*BB 95 : IF(AA*CC.GT.0.) GO TO 3 96 : C 97 : CCCC CC EST NEGATIF ON A 2 SOLUTIONS 1 NEGATIVE 1 POSITIVE 98 : C LA BONNE EST LA POSITIVE 99 : C 100 : SEL=(-B+SQRT(DD))/(2.*AA) 101 : GO TO 333 102 : 3 CONTINUE 103 : C 104 : CCCCCCCC ON A 2 SOLUTIONS POSITIVES OU NEGATIVES 105 : C 106 : IF(DD.GT.0.) GO TO 5 107 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 108 : CCCCCCCCCCCCC LA DROITE S0 STOT NE COUPE PAS LE CONVEXE C 109 : C ON CHERCHE LE POINT DU CONVEXE LE C 110 : CCC PLUS PRES DE STOT ON L APPELLE SIGEL C 111 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 112 : 130 CONTINUE 113 : XX=SELAS/(SSTAR+1.E-20) 113 : CALL AEQBX(SIGEL,STOT,XX,IBOU) 114 : SEL=0. 115 : SPLA=1. 115 : CALL AEQBPC(DSIGP,STOT,SIGEL,UN,-UN,IBOU) 116 : GO TO 555 117 : 5 CONTINUE 118 : C 119 : CCCCCCCCCCCCCCCCCCC ON A 2 SOLUTIONS POSITIVES OU NEGATIVES 120 : C 121 : SEL=(-B+SQRT(DD))/(2.*AA) 122 : GO TO 333 123 : 555 CONTINUE 124 : IF(IMAPLA.NE.11.OR.IN.EQ.1) GO TO 666 125 : IF(IN.EQ.0) GO TO 80 126 : IF(ABS(SELMIN).GE.ABS(SEL)) GO TO 666 127 : SEL=SELMIN 128 : IF(ABS(AAA).LT.TST.AND.ABS(BBB).LT.TST) GO TO 110 129 : IF(ABS(AAA).GT.TST.AND.ABS(CCC).GT.TST.AND.(AAA*CCC).GT.ZERO. 130 : . AND.DDD.LT.ZERO) GO TO 130 131 : IN=1 132 : GO TO 333 133 : 80 SELMIN=SEL 134 : IN=-1 135 : GO TO 90 136 : 666 CONTINUE 137 : RETURN 138 : END 139 :
© Cast3M 2003 - All rights reserved.
Disclaimer