PROGRAM SIGMD
      REAL*8 E1,NU1,P1,D1,C1,R1,E2,NU2,P2,D2,C2,R2,F,FM,D,T,D1M,D2M,
     +P1M,P2M,AE,BE,CE,DE,W,WW,WE,WS,Q1,Q2,Q3,Q4,R3,R4,T3,D3,D4,D3M,
     +D4M,C3,C4,T1B,T1,R,DI,SI,DU,SU,DUL,SUL,DIL,SIL,WSL,WX,RP
      REAL*8 M(3,3),C(3),X(3),FB,FBM
      INTEGER I,J,K
      DATA P1,P2/2.28,0.863/
      DATA D1,D2/525.,1./
      DATA F,D,T/10.,0.1,0.0/
C
C
C CALCUL DE REPARTITION DES CONTRAINTES ET DEFORMATIONS 
C DU SYSTEME : COUCHE MINCE DE SILICE SUR SUBSTRAT DE SILICIUM
C      
C auteur : aime.vareille@wanadoo.fr 
C logiciel libre licence GPL : http://www.gnu.org/copyleft/gpl.html
C      
C E et NU sont les modules dYoung et Coefficient de Poisson Si/SiO2 :
C P1=E1/(1.-NU1) POUR Si (100)
C P2=E2/(1.-NU2) POUR SiO2
C Pour les constantes mécaniques voir par exemple :
C "Intrinsic SiO2 film stress measurements on thermally oxidized Si", J.
C Vac. Sci. Technol. B5(1), pp 15-19, jan/feb 1987, E. Kobeda and E. A.
C Irene.
C
C Le système est considéré de symétrie sphérique :
C Le rayon de courbure moyen du subtrat de Silicium est C1.      
C Le rayon de de courbure moyen de la couche de Silice est C2.      
C Le rayon de la surface neutre du subtrat de Silicium est R1.      
C Le rayon de la surface neutre de la couche de Silice est R2.      
C
C
1     WRITE(6,100)F
100   FORMAT(' Fleche : ',F9.3,' um ? ',$)
      READ(5,*)F
      FM=F*1.E-6
      WRITE(6,101)D
101   FORMAT(' Diametre : ',F9.4,' m ? ',$)
      READ(5,*)D
C      WRITE(6,102)T
102   FORMAT(' Moment flechissant : ',E14.7,' MN ? ',$)
C      READ(5,*)T
      WRITE(6,103)D1
103   FORMAT(' Epaisseur de silicium : ',F9.3,' um ? ',$)
      READ(5,*)D1
      D1M=D1*1.E-6
      WRITE(6,104)D2
104   FORMAT(' Epaisseur de silice : ',F9.3,' um ? ',$)
      READ(5,*)D2
      D2M=D2*1.E-6
      WRITE(6,105)P1
105   FORMAT(' E/(1.-NU) Si (111) 2.28 (100) 1.805: ',F9.5,' E5MPa ? ',
     +$)
      READ(5,*)P1
      WRITE(6,106)P2
106   FORMAT(' E/(1.-NU) SiO2 0.863 MPa : ',F9.5,' E5MPa ? ',$)
      READ(5,*)P2
      P1M=P1*1.E5
      P2M=P2*1.E5
C
C RAYON DE COURBURE EN SURFACE R
C
      R=(FM**2+D**2/4.)/2./FM
      Q1=P1M*D1M
      Q2=P2M*D2M
      C1=R-D2M-D1M/2.
      C2=R-D2M/2.
C
C
C LIGNE NEUTRE DU SILICIUM (P1) SUR LE RAYON R1
C MOMENT DE FLEXION T1
C LIGNE NEUTRE DE LA SILICE (P2) SUR LE RAYON R2
C
C
      T1B=Q1*(C2-C1)+Q2*D2M*D2M/12./C2*(1.+Q1/Q2)
      T1=T+T1B
      R1=Q1*((C2-C1)*C1+D2M*D2M*C1/12./C2-D1M*D1M/12.)/T1
      R2=C2/(1.-Q1/Q2*(C1/R1-1.))
C
C
C CONTRAINTE DANS LE SILICIUM A L'INTERFACE SI A LA SURFACE SU
C DEFORMATION DANS LE SILICIUM A L'INTERFACE DI A LA SURFACE DU
C CONTRAINTE DANS LA SILICE A L'INTERFACE SIL A LA SURFACE SUL
C DEFORMATION DANS LA SILICE A L'INTERFACE DIL A LA SURFACE DUL
C
C
      DI=(R-D2M-R1)/R1
      SI=P1M*DI
      DU=(R-D2M-D1M-R1)/R1
      SU=P1M*DU
      DUL=(R-R2)/R2
      SUL=P2M*DUL
      DIL=(R-D2M-R2)/R2
      SIL=P2M*DIL
C
C
C ENERGIE ELASTIQUE DU SILICIUM WS
C ENERGIE ELASTIQUE DE LA SILICE WSL
C
C
      WS=Q1/R1/R1*((C1-R1)**2+D1M*D1M/12.)*1.E6
      WSL=Q2/R2/R2*((C2-R2)**2+D2M*D2M/12.)*1.E6
      WRITE(6,222)F,D,D1,D2
222   FORMAT(' Fleche : ',F9.3,' um ; diametre : ',F9.4,' m'/
     +' Epaisseur de silicium : ',F9.3,' um, de silice : ',F9.3,' um')
      WRITE(6,200)R,T1,R1-R,R2-R,SI,SU,SIL,SUL,DI,DU,DIL,DUL,WS,WSL
200   FORMAT(' Rayon de courbure : ',F9.3,'m ; moment de flexion : ',
     +E14.7,'MN'/' Lignes neutre du silicium : ',E14.7,
     +'m, de la silice : ',E14.7,'m'/
     +' Contrainte dans le silicium : interface ',E14.7,' surface : ',
     +E14.7,'MPa'/' Contrainte dans la silice : interface ',E14.7,
     +' surface : ',E14.7,'MPa'/
     +' Deformation silicium : interface ',E14.7,' surface : ',
     +E14.7/' Deformation silice : interface ',E14.7,
     +' surface : ',E14.7/' Energie elastique : silicium ',E14.7,
     +', silice ',E14.7,' J.m-2',/)
      K=0
      WRITE(6,224)K
224   FORMAT(' Continue ou recommence (0/1) : ',I2,' ? ',$)
      READ(5,*)K
       IF(K.EQ.1) GOTO 1
C
C
C UN AMINCICEMENT RETIRE AU SYSTEME L'ENERGIE ELASTIQUE CORRESPONDANTE
C LA RESOLUTION UTILISE LA METHODE DE NEWTON :
C CALCUL DE L'ENERGIE ELASTIQUE RESTANTE
C           PUIS CALCUL DU JACOBIEN
C           JUSQU'A INVARIABILITE DES VALEURS DE C1, C2, R1 ET R2
C
C
      D3=D1
      D4=D2
2     WRITE(6,110)D3
110   FORMAT(' Epaisseur de silicium : ',F9.3,' um ? ',$)
      READ(5,*)D3
      D3M=D3*1.E-6
      WRITE(6,111)D4
111   FORMAT(' Epaisseur de silice : ',F9.3,' um ? ',$)
      READ(5,*)D4
      TE=T*1.E6
      WRITE(6,112)TE
112   FORMAT(' Moment flechissant : ',E14.7,' N ? ',$)
      READ(5,*)TE
      T=TE*1.E-6
      D4M=D4*1.E-6
      Q3=P1M*D3M
      Q4=P2M*D4M
      W=Q3/R1**2*((C1-R1)**2+D3M**2/12.)+
     +Q4/R2**2*((C2-R2)**2+D4M**2/12.)
      C3=C1
      C4=C2
      R3=R1
      R4=R2
      K=0
C
C CALCUL DU JACOBIEN LITTERAL !!!
C
3     M(1,1)=Q3*R4+Q4*R3
      M(2,1)=2.*Q3*R4*(C3-R3)+Q4*R3*(2.*C3-R3-R4+D3M+D4M)
      M(3,1)=2.*Q3*R4**2*(C3-R3)+2.*Q4*R3**2*(C3-R4+(D3M+D4M)/2.)
      M(1,2)=-1.*Q3*R4+Q4*(C3-R4+(D3M+D4M)/2.)
      M(2,2)=2.*Q3*R4*(R3-C3)-Q4*R3*(C3-R4+(D3M+D4M)/2.)+
     +Q4*((C3-R4+(D3M+D4M)/2.)*(C3-R3+(D3M+D4M)/2.)+D4M**2/12.)+
     +T*R4 
      M(3,2)=2.*Q3*R4**2*(R3-C3)+2.*Q4*R3*((C3-R4+(D3M+D4M)/2.)+
     +D4M**2/12.)-2.*W*R3*R4**2
      M(1,3)=Q3*(C3-R3)-Q4*R3
      M(2,3)=Q3*((C3-R3)**2+D3M**2/12.)-Q4*R3*(C3-R3+(D3M+D4M)/2.)+
     +T*R3 
      M(3,3)=2.*Q3*R4*((C3-R3)**2+D3M**2/12.)-
     +2.*Q4*R3**2*(C3-R4+(D3M+D4M)/2.)-2.*W*R3**2*R4
      K=K+1
      C(1)=M(1,1)*C3+M(1,2)*R3+M(1,3)*R4-Q3*R4*(C3-R3)-
     +Q4*R3*(C3-R4+(D3M+D4M)/2.)
      C(2)=M(2,1)*C3+M(2,2)*R3+M(2,3)*R4-Q3*R4*((C3-R3)**2+D3M**2/12.)-
     +Q4*R3*((C3-R4+(D3M+D4M)/2.)*(C3-R3+(D3M+D4M)/2.)+D4M**2/12.)-
     +T*R3*R4 
      C(3)=M(3,1)*C3+M(3,2)*R3+M(3,3)*R4-
     +Q3*R4**2*((C3-R3)**2+D3M**2/12.)-
     +Q4*R3**2*((C3-R4+(D3M+D4M)/2.)**2+D4M**2/12.)+W*R3**2*R4**2
      CALL GAUS3(M,C,X)
      C3=X(1)
      R3=X(2)
      R4=X(3)
      C4=C3+(D3M+D4M)/2.
      WX=Q3/R3**2*((C3-R3)**2+D3M**2/12.)+
     +Q4/R4**2*((C4-R4)**2+D4M**2/12.)
C      WRITE(6,201)C1,C2,R1,R2,C3,C4,R3,R4,W,WX
201   FORMAT(' C1,C2,R1,R2 :',4(1X,E14.7)/
     +' C3,C4,R3,R4 :',4(1X,E14.7)/
     +' W,WX :',2(1X,E14.7)/)
C
C ARRET DES ITERATIONS LORSQUE L'ECART D'ENERGIE ATTEINT MOINS DE 1/10000 OU
C                      LORSQUE LE NOMBRE DES ITERATIONS EST SUPERIEUR A 10000
C
       IF((ABS(W-WX).LT.W/10000.).OR.(K.EQ.10000)) THEN
      WRITE(6,201)C1,C2,R1,R2,C3,C4,R3,R4,W,WX
C
C
C CONTRAINTE DANS LE SILICIUM A L'INTERFACE SI A LA SURFACE SU
C DEFORMATION DANS LE SILICIUM A L'INTERFACE DI A LA SURFACE DU
C CONTRAINTE DANS LA SILICE A L'INTERFACE SIL A LA SURFACE SUL
C DEFORMATION DANS LA SILICE A L'INTERFACE DIL A LA SURFACE DUL
C
C
      DI=(C4-D4M/2.-R3)/R3
      SI=P1M*DI
      DU=(C4-D4M/2.-D3M-R3)/R3
      SU=P1M*DU
      DUL=(C4+D4M/2.-R4)/R4
      SUL=P2M*DUL
      DIL=(C4-D4M/2.-R4)/R4
      SIL=P2M*DIL
C
C
C ENERGIE ELASTIQUE DU SILICIUM WS
C ENERGIE ELASTIQUE DE LA SILICE WSL
C
C
      WS=Q3/R3**2*((C3-R3)**2+D3M*D3M/12.)*1.E6
      WSL=Q4/R4**2*((C4-R4)**2+D4M*D4M/12.)*1.E6
      RP=C4+D4M/2.
      T1B=Q3*(C4-C3)+Q4*D4M*D4M/12./C4*(1.+Q3/Q4)
      T1=T+T1B
      AE=1.0
      BE=-2.*RP
      CE=D**2/4.
      DE=BE**2-4.*AE*CE
      FB=1.E6*((-BE-SIGN(SQRT(DE),RP))/2./AE)
      WRITE(6,223)FB,D,D3,D4,K
223   FORMAT(' Fleche : ',F9.3,' um ; diametre : ',F9.4,' m'/
     +' Epaisseur de silicium : ',F9.3,' um, de silice : ',F9.3,' um'/
     +1X,I5,' regressions')
      WRITE(6,200)RP,T1,R3-RP,R4-RP,SI,SU,SIL,SUL,DI,DU,DIL,DUL,WS,WSL
C200   FORMAT(' Rayon de courbure : ',F9.3,'m ; moment de flexion : ',
C     +E14.7,'MN'/' Lignes neutre du silicium : ',E14.7,
C     +'m, de la silice : ',E14.7,'m'/
C     +' Contrainte dans le silicium : interface ',F9.4,' surface : ',
C     +F9.4,'MPa'/' Contrainte dans la silice : interface ',F9.4,
C     +' surface : 'F10.4,'MPa'/
C     +' Deformation silicium : interface ',E14.7,' surface : ',
C     +E14.7/' Deformation silice : interface ',E14.7,
C     +' surface : 'E14.7/' Energie elastique : silicium ',E14.7,
C     +', silice ',E14.7,' J.m-2',/)
      K=0
      WRITE(6,224)K
C224   FORMAT(' Continue ou recommence (0/1) : ',I2,' ? ')
      READ(5,*)K
       IF(K.EQ.1) GOTO 1
       GOTO 2
       ENDIF
      GOTO 3
      END
	SUBROUTINE GAUS3(M3,C3,X)
	REAL*8 X(3),
     +M3(3,3),M2(2,2),M,C3(3),C2(2),C
	INTEGER I,J,K,L,A
C
C METHODE DE GAUS DE RESOLUTION DES SYSTEMES LINEAIRES
C
C
C 1) REMPLISSAGE DES TABLEAUX M1-2
C 2) RESOLUTION
C
C	WRITE(6,203)((M3(I,J),J=1,3),I=1,3)
203	FORMAT(3(3(1X,E14.7)/))
C	ACCEPT *,A
	  DO 7 I=1,2
	  C2(I)=C3(I)*M3(3,3)-C3(3)*M3(I,3)
	    DO 8 J=1,2
	    M2(I,J)=M3(I,J)*M3(3,3)-M3(3,J)*M3(I,3)
8	    CONTINUE
7	  CONTINUE
C	WRITE(6,202)((M2(I,J),J=1,2),I=1,2)
202	FORMAT(2(2(1X,E14.7)/))
C	ACCEPT *,A
	C=C2(1)*M2(2,2)-C2(2)*M2(1,2)
        M=M2(1,1)*M2(2,2)-M2(2,1)*M2(1,2)
C RESOLUTION
	  IF(M.EQ.0) GO TO 91
	X(1)=C/M
	  IF(M2(2,2).NE.0) THEN
	  X(2)=(C2(2)-M2(2,1)*X(1))/M2(2,2)
	  ELSE
	    IF(M2(1,2).EQ.0) GO TO 92
	  X(2)=(C2(1)-M2(1,1)*X(1))/M2(1,2)
	  ENDIF
	  IF(M3(3,3).NE.0) THEN
	  X(3)=(C3(3)-M3(3,1)*X(1)-M3(3,2)*X(2))/M3(3,3)
	  ELSE
	    IF(M3(2,3).NE.0) THEN
	    X(3)=(C3(2)-M3(2,1)*X(1)-M3(2,2)*X(2))/M3(2,3)
	    ELSE
	      IF(M3(1,3).EQ.0) GO TO 93
	    X(3)=(C3(1)-M3(1,1)*X(1)-M3(1,2)*X(2))/M3(1,3)
	    ENDIF
	  ENDIF
	RETURN
91	WRITE(6,101)
101	FORMAT(' GAUS3 ne peut pas resoudre.')
	RETURN
92	WRITE(6,102)X(1)
102	FORMAT(' GAUS3 ne peut pas resoudre: X(1) = ',E14.7)
	RETURN
93	WRITE(6,103)X(1),X(2)
103	FORMAT(' GAUS3 ne peut pas resoudre:'/' X(1-2) = ',2(1X,E14.7))
	RETURN
	END