PROGRAM cijrln C----------------------------------------------------- C EULER AVERAGE OF C(i,j)*R(l,niu) FOR l=0,2,4 C necesary to calculate mean stress tensor in the WSODF c representation of Voigt ground state. c Input file: 'CSC.TXT': SINGLE CRYSTAL STIFFNESS CONSTANTS C Output file: 'CRIJLN.TXT' C----------------------------------------------------- COMMON /bij/i,j common/lniu/l,niu common /cmono/c(6,6) DIMENSION CRLN(6,6) EXTERNAL CR OPEN(1,FILE='CSC.TXT') DO 1 I=1,6 READ(1,*)(C(I,J),J=1,6) 1 CONTINUE CLOSE(1) OPEN(1,FILE='CRIJLN.TXT') 200 FORMAT(1X,6(F9.6,2X)) PRINT 100 100 FORMAT(/' ENTER NS >= 21 - odd - for SIMPSON'/) READ(5,*)NS DO 3 L=0,4,2 NL=(2*L+1)**2 DO 3 niu=1,NL call gmnv DO 5 I=1,6 DO 5 J=i,6 CALL INTEUL(NS,CR,CRLN(I,J)) 5 CONTINUE do 500 j=1,5 do 500 i=j+1,6 500 crln(i,j)=crln(j,i) PRINT 95,L,Niu WRITE(1,95)L,Niu 95 FORMAT(/10X,' L= ',I2,' NIU = ',I2/) DO 4 J=1,6 PRINT 200,(CRLN(I,J),I=1,6) WRITE(1,200)(CRLN(I,J),I=1,6) 4 CONTINUE 3 CONTINUE CLOSE(1) END C --------------------------------------------------- subroutine gmnv c Correspondence between g(i,l,niu) and c (alpha,beta,gamma,delta)(i,l,m,n) c iabgd = 1, 2, 3, 4 common /lniu/l,niu common /mni/m,n,iabgd l21=2*l+1 1 if(mod(niu,l21).eq.0)then nius=niu/l21 niuc=l21 else nius=niu/l21+1 niuc=niu-(nius-1)*l21 endif if(nius.eq.1.or.mod(nius,2).eq.0)go to 1000 c gama & delta if(niuc.eq.1.or.mod(niuc,2).eq.0)then iabgd=3 else iabgd=4 endif go to 3000 c alpha,beta 1000 if(niuc.eq.1.or.mod(niuc,2).eq.0)then iabgd=1 else iabgd=2 endif 3000 m=niuc/2 n=nius/2 return end c ----------------------------------------------------- FUNCTION CR(F1,F2,X) COMMON /bij/i,j common /lniu/l,niu CR=CIJ(I,J,F1,X,F2)*rlniu(F1,X,F2) RETURN END C ----------------------------------------------------------- FUNCTION CIJ(I,J,F1,X,F2) C SINGLE CRYSTAL ELASTIC CONSTANTS IN SAMPLE REFERENCE SYSTEM C formula (74a) din carte common /cmono/c(6,6) dimension ro(6) data ro/1.,1.,1.,2.,2.,2./ CIJ=0. DO 1 L=1,6 W=0. DO 2 K=1,6 W=W+C(L,K)*RO(K)*Q(K,J,F1,X,F2) 2 CONTINUE CIJ=CIJ+P(I,L,F1,X,F2)*W 1 CONTINUE CIJ=CIJ/RO(J) RETURN END C ------------------------------------------ FUNCTION A(I,J,F1,X,F2) C THE EULER'S MATRIX A(3,3) GO TO (1,2,3),I 1 GO TO (11,12,13),J 11 A=COS(F1)*COS(F2)-SIN(F1)*SIN(F2)*X RETURN 12 A=SIN(F1)*COS(F2)+COS(F1)*SIN(F2)*X RETURN 13 A=SIN(F2)*SQRT(1.-X*X) RETURN 2 GO TO (21,22,23),J 21 A=-COS(F1)*SIN(F2)-SIN(F1)*COS(F2)*X RETURN 22 A=-SIN(F1)*SIN(F2)+COS(F1)*COS(F2)*X RETURN 23 A=COS(F2)*SQRT(1.-X*X) RETURN 3 GO TO (31,32,33),J 31 A=SIN(F1)*SQRT(1.-X*X) RETURN 32 A=-COS(F1)*SQRT(1.-X*X) RETURN 33 A=X RETURN END C ---------------------------------------- FUNCTION P(I,J,F1,X,F2) C THE MATRIX P(6,6) GO TO (1,2,3,4,5,6),I 1 GO TO (11,12,13,14,15,16),J 11 P=A(1,1,F1,X,F2)**2 RETURN 12 P=A(2,1,F1,X,F2)**2 RETURN 13 P=A(3,1,F1,X,F2)**2 RETURN 14 P=2*A(2,1,F1,X,F2)*A(3,1,F1,X,F2) RETURN 15 P=2*A(1,1,F1,X,F2)*A(3,1,F1,X,F2) RETURN 16 P=2*A(1,1,F1,X,F2)*A(2,1,F1,X,F2) RETURN 2 GO TO (21,22,23,24,25,26),J 21 P=A(1,2,F1,X,F2)**2 RETURN 22 P=A(2,2,F1,X,F2)**2 RETURN 23 P=A(3,2,F1,X,F2)**2 RETURN 24 P=2*A(2,2,F1,X,F2)*A(3,2,F1,X,F2) RETURN 25 P=2*A(1,2,F1,X,F2)*A(3,2,F1,X,F2) RETURN 26 P=2*A(1,2,F1,X,F2)*A(2,2,F1,X,F2) RETURN 3 GO TO (31,32,33,34,35,36),J 31 P=A(1,3,F1,X,F2)**2 RETURN 32 P=A(2,3,F1,X,F2)**2 RETURN 33 P=A(3,3,F1,X,F2)**2 RETURN 34 P=2*A(2,3,F1,X,F2)*A(3,3,F1,X,F2) RETURN 35 P=2*A(1,3,F1,X,F2)*A(3,3,F1,X,F2) RETURN 36 P=2*A(1,3,F1,X,F2)*A(2,3,F1,X,F2) RETURN 4 GO TO (41,42,43,44,45,46),J 41 P=A(1,2,F1,X,F2)*A(1,3,F1,X,F2) RETURN 42 P=A(2,2,F1,X,F2)*A(2,3,F1,X,F2) RETURN 43 P=A(3,2,F1,X,F2)*A(3,3,F1,X,F2) RETURN 44 P=A(2,2,F1,X,F2)*A(3,3,F1,X,F2)+A(3,2,F1,X,F2)*A(2,3,F1,X,F2) RETURN 45 P=A(1,2,F1,X,F2)*A(3,3,F1,X,F2)+A(3,2,F1,X,F2)*A(1,3,F1,X,F2) RETURN 46 P=A(1,2,F1,X,F2)*A(2,3,F1,X,F2)+A(2,2,F1,X,F2)*A(1,3,F1,X,F2) RETURN 5 GO TO (51,52,53,54,55,56),J 51 P=A(1,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 52 P=A(2,1,F1,X,F2)*A(2,3,F1,X,F2) RETURN 53 P=A(3,1,F1,X,F2)*A(3,3,F1,X,F2) RETURN 54 P=A(2,1,F1,X,F2)*A(3,3,F1,X,F2)+A(3,1,F1,X,F2)*A(2,3,F1,X,F2) RETURN 55 P=A(1,1,F1,X,F2)*A(3,3,F1,X,F2)+A(3,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 56 P=A(1,1,F1,X,F2)*A(2,3,F1,X,F2)+A(2,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 6 GO TO (61,62,63,64,65,66),J 61 P=A(1,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN 62 P=A(2,1,F1,X,F2)*A(2,2,F1,X,F2) RETURN 63 P=A(3,1,F1,X,F2)*A(3,2,F1,X,F2) RETURN 64 P=A(2,1,F1,X,F2)*A(3,2,F1,X,F2)+A(3,1,F1,X,F2)*A(2,2,F1,X,F2) RETURN 65 P=A(1,1,F1,X,F2)*A(3,2,F1,X,F2)+A(3,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN 66 P=A(1,1,F1,X,F2)*A(2,2,F1,X,F2)+A(2,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN END C ---------------------------------------------------------------------- FUNCTION Q(I,J,F1,X,F2) C THE MATRIX Q(6,6) GO TO (1,2,3,4,5,6),I 1 GO TO (11,12,13,14,15,16),J 11 Q=A(1,1,F1,X,F2)**2 RETURN 12 Q=A(1,2,F1,X,F2)**2 RETURN 13 Q=A(1,3,F1,X,F2)**2 RETURN 14 Q=2*A(1,2,F1,X,F2)*A(1,3,F1,X,F2) RETURN 15 Q=2*A(1,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 16 Q=2*A(1,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN 2 GO TO (21,22,23,24,25,26),J 21 Q=A(2,1,F1,X,F2)**2 RETURN 22 Q=A(2,2,F1,X,F2)**2 RETURN 23 Q=A(2,3,F1,X,F2)**2 RETURN 24 Q=2*A(2,2,F1,X,F2)*A(2,3,F1,X,F2) RETURN 25 Q=2*A(2,1,F1,X,F2)*A(2,3,F1,X,F2) RETURN 26 Q=2*A(2,1,F1,X,F2)*A(2,2,F1,X,F2) RETURN 3 GO TO (31,32,33,34,35,36),J 31 Q=A(3,1,F1,X,F2)**2 RETURN 32 Q=A(3,2,F1,X,F2)**2 RETURN 33 Q=A(3,3,F1,X,F2)**2 RETURN 34 Q=2*A(3,2,F1,X,F2)*A(3,3,F1,X,F2) RETURN 35 Q=2*A(3,1,F1,X,F2)*A(3,3,F1,X,F2) RETURN 36 Q=2*A(3,1,F1,X,F2)*A(3,2,F1,X,F2) RETURN 4 GO TO (41,42,43,44,45,46),J 41 Q=A(2,1,F1,X,F2)*A(3,1,F1,X,F2) RETURN 42 Q=A(2,2,F1,X,F2)*A(3,2,F1,X,F2) RETURN 43 Q=A(2,3,F1,X,F2)*A(3,3,F1,X,F2) RETURN 44 Q=A(2,2,F1,X,F2)*A(3,3,F1,X,F2)+A(3,2,F1,X,F2)*A(2,3,F1,X,F2) RETURN 45 Q=A(2,1,F1,X,F2)*A(3,3,F1,X,F2)+A(3,1,F1,X,F2)*A(2,3,F1,X,F2) RETURN 46 Q=A(2,1,F1,X,F2)*A(3,2,F1,X,F2)+A(2,2,F1,X,F2)*A(3,1,F1,X,F2) RETURN 5 GO TO (51,52,53,54,55,56),J 51 Q=A(1,1,F1,X,F2)*A(3,1,F1,X,F2) RETURN 52 Q=A(1,2,F1,X,F2)*A(3,2,F1,X,F2) RETURN 53 Q=A(1,3,F1,X,F2)*A(3,3,F1,X,F2) RETURN 54 Q=A(1,2,F1,X,F2)*A(3,3,F1,X,F2)+A(3,2,F1,X,F2)*A(1,3,F1,X,F2) RETURN 55 Q=A(1,1,F1,X,F2)*A(3,3,F1,X,F2)+A(3,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 56 Q=A(1,1,F1,X,F2)*A(3,2,F1,X,F2)+A(3,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN 6 GO TO (61,62,63,64,65,66),J 61 Q=A(1,1,F1,X,F2)*A(2,1,F1,X,F2) RETURN 62 Q=A(1,2,F1,X,F2)*A(2,2,F1,X,F2) RETURN 63 Q=A(1,3,F1,X,F2)*A(2,3,F1,X,F2) RETURN 64 Q=A(1,2,F1,X,F2)*A(2,3,F1,X,F2)+A(1,3,F1,X,F2)*A(2,2,F1,X,F2) RETURN 65 Q=A(1,1,F1,X,F2)*A(2,3,F1,X,F2)+A(2,1,F1,X,F2)*A(1,3,F1,X,F2) RETURN 66 Q=A(1,1,F1,X,F2)*A(2,2,F1,X,F2)+A(2,1,F1,X,F2)*A(1,2,F1,X,F2) RETURN END C -------------------------------------------------------------------- function rlniu(f1,x,f2) C Functions R(l,niu) common /lniu/l,niu common /mni/m,n,iabgd mn=m+n m0=(-1)**(l+n) go to (1,2,3,4),iabgd 1 if(mod(mn,2).eq.0)then rlniu=(cos(m*f2+n*f1)*qlmn(l,m,n,x)+m0* $ cos(m*f2-n*f1)*qlmn(l,m,n,-x))/2. else rlniu=(sin(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ sin(m*f2-n*f1)*qlmn(l,m,n,-x))/2. endif return 2 if(mod(mn,2).eq.0)then rlniu=(-sin(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ sin(m*f2-n*f1)*qlmn(l,m,n,-x))/2. else rlniu=(cos(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ cos(m*f2-n*f1)*qlmn(l,m,n,-x))/2. endif return 3 if(mod(mn,2).eq.0)then rlniu=(sin(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ sin(m*f2-n*f1)*qlmn(l,m,n,-x))/2. else rlniu=(-cos(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ cos(m*f2-n*f1)*qlmn(l,m,n,-x))/2. endif return 4 if(mod(mn,2).eq.0)then rlniu=(cos(m*f2+n*f1)*qlmn(l,m,n,x)-m0* $ cos(m*f2-n*f1)*qlmn(l,m,n,-x))/2. else rlniu=(sin(m*f2+n*f1)*qlmn(l,m,n,x)+m0* $ sin(m*f2-n*f1)*qlmn(l,m,n,-x))/2. endif return end C ------------------------------------------------------ FUNCTION qlmn(l,m,n,z) c Real Legendre functions: c Qlmn=Plmn for m+n=even and Qlmn=i*Plmn for m+n=odd if(l.gt.0)go to 1000 c l=0 ======= Qlmn=1. return 1000 mn=m+n go to (1,2,3,4),l c l=1 ======= 1 if(mn.gt.0)go to 1001 c Q100 Qlmn=z return c Q110=Q101 & Q111 1001 if(mn.eq.1)then qlmn=0.7071067811865475*sqrt(1.-z**2) else qlmn=(1.+z)/2. endif return c l=2 ======= 2 if(mn.gt.0)go to 1002 c Q200 qlmn=(3*Z**2-1)/2 return 1002 if(m.eq.2.or.n.eq.2)go to 1003 c Q210=Q201 & Q211 if(mn.eq.1)then qlmn=1.224744871392*Z*SQRT(1.-Z**2) else qlmn=(2*Z**2+Z-1.)/2 endif return c Q220=Q202 & Q221=Q212 & Q222 1003 if(mn.eq.2)then qlmn=-0.6123724356958*(1.-Z**2) elseif(mn.eq.3)then qlmn=(1.+Z)*SQRT(1.-Z**2)/2 else qlmn=(1.+Z)**2/4 endif return c l=3 ======= 3 if(mn.gt.0)go to 1004 c Q300 qlmn=z*(5*z**2-3.)/2 return 1004 if(m.gt.1.or.n.gt.1)go to 1005 c Q310=Q301 & Q311 if(mn.eq.1)then qlmn=0.4330127018922193*sqrt(1.-z**2)*(5*z**2-1.) else qlmn=(1.+z)*(15*z**2-10*z-1.)/8 endif return 1005 if(m.eq.3.or.n.eq.3)go to 1006 c Q320=Q302 & Q321=Q312 & Q322 if(mn.eq.2)then qlmn=-1.369306393762915*Z*(1.-z**2) elseif(mn.eq.3)then qlmn=0.3952847075210474*sqrt(1.-z**2)*(1.+z)*(3*z-1.) else qlmn=(1.+z)**2*(3*z-2)/4 endif return c Q330=Q303 & Q331=Q313 & Q332=Q323 & Q333 1006 if(mn.eq.3)then qlmn=-0.5590169943749474*(1.-z**2)**1.5 elseif(mn.eq.4)then qlmn=-0.4841229182759271*(1.-z)*(1.+z)**2 elseif(mn.eq.5)then qlmn=0.3061862178478973*sqrt(1.-z**2)*(1.+z)**2 else qlmn=(1.+z)**3/8 endif return c l=4 ======== 4 if(mn.gt.0)go to 1007 c Q400 qlmn=(35*Z**4-30*Z**2+3)/8. return 1007 if(m.gt.1.or.n.gt.1)go to 1008 c Q410=Q401 & Q411 if(mn.eq.1)then qlmn=0.5590169943749*Z*SQRT(1.-Z**2)*(7*Z**2-3.) else qlmn=(1.+Z)*(28*Z**3-21*Z**2-6*Z+3.)/8. endif return 1008 if(m.gt.2.or.n.gt.2)go to 1009 c Q420=q402 & Q421=Q412 & Q422 If(mn.eq.2)then qlmn=-0.395284707521*(1.-Z**2)*(7*Z**2-1.) elseif(mn.eq.3)then qlmn=0.1767766952966*SQRT(1.-Z**2)*(1.+Z)*(14*Z**2-7*Z-1.) else qlmn=(1.+Z)**2*(7*Z**2-7*Z+1.)/4. endif return 1009 if(m.eq.4.or.n.eq.4)go to 1010 c Q430=Q403 & Q431=Q413 & Q432=Q423 & Q433 if(mn.eq.3)then Qlmn=-1.479019945775*Z*(1.-Z**2)**1.5 elseif(mn.eq.4)then Qlmn=-0.330718913883*(1.-Z**2)*(4*Z**2+3*Z-1.) elseif(mn.eq.5)then Qlmn=0.4677071733467*SQRT(1.-Z**2)*(1.+Z)**2*(2*Z-1.) else Qlmn=(1.+Z)**3*(4*Z-3.)/8. endif return c Q440=Q404 & Q441=Q414 & q442=Q424 & Q443=Q434 & Q444 1010 if(mn.eq.4)then qlmn=0.5229125165838*(1.-Z**2)**2 elseif(mn.eq.5)then Qlmn=-0.4677071733467*(1.-Z**2)**1.5*(1.+Z) elseif(mn.eq.6)then Qlmn=-0.330718913883*(1.-Z**2)*(1.+Z)**2 elseif(mn.eq.7)then Qlmn=0.1767766952966*SQRT(1.-Z**2)*(1.+Z)**3 else qlmn=(1.+Z)**4/16. endif return end c --------------------------------------------------------- SUBROUTINE INTEUL(N,FUN,REZ) C INTEGRAL OVER EULER ANGLES. C OVER F1,F2 (0,2PI) INTEGRATED BY SIMPSON N(ODD) NODES C OVER X=COS(F0) (-1,1) INTEGRATED BY GAUSS 10 NODES REAL*8 H DIMENSION CS(101),CG(5),XG(5) DATA CG/.0666713,.1494513,.2190864,.2692667,.2955242/ DATA XG/.9739065,.8650634,.6794096,.4333954,.1488743/ h=8*datan(1.d0)/(N-1) M=(N-1)/2 CS(1)=1. CS(N)=1. DO 1 K=1,M 1 CS(2*K)=4. DO 2 K=1,M-1 2 CS(2*K+1)=2. REZ=0. DO 3 I=1,N F1=(I-1)*H RJ=0. DO 4 J=1,N F2=(J-1)*H RK=0. DO 5 K=1,5 5 RK=RK+CG(K)*(FUN(F1,F2,XG(K))+FUN(F1,F2,-XG(K))) 4 RJ=RJ+CS(J)*RK 3 REZ=REZ+CS(I)*RJ REZ=REZ/18/(N-1)**2 RETURN END