C
C
SUBROUTINE CANORA(XP,YP,X1,Y1,X2,Y2,C,MU1,MU2,QP)
C
C La rutina CANORA Calcula la matriz Cij del termino libre
C analiticamente segun la formulacion de Vlado, para materiales
c ortótropos fuera de ejes de ortotropía y anisótropos.
C
IMPLICIT NONE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C VARIABLES EXTERNAS (ENTRADA)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C MU1,MU2 : Solucion del polinomio caracteristico
COMPLEX*16 MU1,MU2
COMPLEX*16 QP(2,2,2)
C XP,YP : Coordenadas del punto de colocacion
REAL*8 XP,YP
C X1,Y1,X2,Y2 : Coordenadas de los nodos anterior y posterior
REAL*8 X1,Y1,X2,Y2
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C VARIABLES EXTERNAS (SALIDA)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C C(,) : Termino libre
REAL*8 C(2,2)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCC
C VARIABLES INTERNAS
CCCCCCCCCCCCCCCCCCCCCCCCCCC
C
REAL*8 R(2,2),RL1(2,2),RL2(2,2)
COMPLEX*16 IC,Z,Z1(2),Z2(2),PRUEBA
REAL*8 XB,XA,YB,YA,X,Y,PI,EL1,EL2
REAL*8 ARG0(2),RMOD,DIS,CMOD
INTEGER I,J,L
C
C LA FUNCION CMOD() CALCULA EL MODULO DE UN COMPLEJO
C LA FUNCION DIS() CALCULA LA DISTANCIA ENTRE DOS PUNTOS
C LA FUNCION RMOD() CALCULA EL MODULO DE UN VECTOR
CMOD(Z)=DSQRT(DREAL(Z)*DREAL(Z)+DIMAG(Z)*DIMAG(Z))
DIS(XA,YA,XB,YB)=DSQRT((XB-XA)*(XB-XA)+(YB-YA)*(YB-YA))
RMOD(X,Y)=DSQRT(X*X+Y*Y)
C
CCCCCCCCCCCCCCC
C RUTINA
CCCCCCCCCCCCCCC
C
PI=4.D00*DATAN(1.D00)
IC=(0.D00,1.D00)
EL1=DIS(XP,YP,X1,Y1)
EL2=DIS(XP,YP,X2,Y2)
C
C R(I,J) ES LA COORDENADA J DEL VECTOR UNITARIO TANGENTE AL
C CONTORNO (DISCRETIZADO) POR EL LADO I
C
R(1,1)=(X1-XP)/EL1
R(1,2)=(Y1-YP)/EL1
R(2,1)=(X2-XP)/EL2
R(2,2)=(Y2-YP)/EL2
C
C RL1(L,J) ES LA COORDENADA J DEL VECTOR R1 ASOCIADO A LA RAIZ
C L DEL POLINOMIO CARACTERISTICO
C
RL1(1,1)=R(1,1)+DREAL(MU1)*R(1,2)
RL1(1,2)=DIMAG(MU1)*R(1,2)
C IF (.NOT.(MU1.EQ.MU2)) THEN
RL1(2,1)=R(1,1)+DREAL(MU2)*R(1,2)
RL1(2,2)=DIMAG(MU2)*R(1,2)
C END IF
C
RL2(1,1)=R(2,1)+DREAL(MU1)*R(2,2)
RL2(1,2)=DIMAG(MU1)*R(2,2)
C IF (.NOT.(MU1.EQ.MU2)) THEN
RL2(2,1)=R(2,1)+DREAL(MU2)*R(2,2)
RL2(2,2)=DIMAG(MU2)*R(2,2)
C
C
C Z1=DCMPLX(RL1(1,1),RL1(1,2))
C Z2=DCMPLX(RL2(1,1),RL2(1,2))
C
C ARG0=DATAN2(DIMAG(Z1/Z2),DREAL(Z1/Z2))
C IF (ARG0.LT.0.D00) THEN
C ARG0=2.*PI+ARG0
C END IF
C
C DO I=1,2
C DO J=1,2
C C(I,J)=DREAL(-T1(I,J)*(DLOG(CMOD(Z1)/CMOD(Z2))+IC*
C * ARG0)+T2(I,J)*(CONJG(Z2)/Z2-CONJG(Z1)/Z1))
C END DO
C END DO
C END IF
C END IF
C RETURN
C
C ESTOS SON LOS TRANSFORMADOS COMPLEJOS PARA LA RAIZ MU1
C
Z1(1)=DCMPLX(RL1(1,1),RL1(1,2))
Z2(1)=DCMPLX(RL2(1,1),RL2(1,2))
C
C ESTOS SON LOS TRANSFORMADOS COMPLEJOS PARA LA RAIZ MU2
C
Z1(2)=DCMPLX(RL1(2,1),RL1(2,2))
Z2(2)=DCMPLX(RL2(2,1),RL2(2,2))
C
C if((dabs(dimag(z1/z2)))).lt.1.d-14)then
C arg0=pi
C else
C arg0=pi*(1.d0-dsign(1.d0,(-rl2(1,2)*rl1(1,1)+rl2(1,1)*
C @rl1(1,2))))+datan2((-rl2(1,2)*rl1(1,1)+rl2(1,1)*rl1(1,2)),
C @ (rl2(1,1)*rl1(1,1)+rl2(1,2)*rl1(1,2)))
C end if
C
PRUEBA=LOG(Z1(1)/Z2(1))
C IF(DIMAG(Z1(1)/Z2(1)).LT.0.d00) THEN
C ARG0(1)=2.D00*PI+DATAN2(DIMAG(Z1(1)/Z2(1)),DREAL(Z1(1)/Z2(1)))
C ELSE
C ARG0(1)=DATAN2(DIMAG(Z1(1)/Z2(1)),DREAL(Z1(1)/Z2(1)))
C END IF
IF(DIMAG (LOG(Z1(1)/Z2(1))).LT.0.d00) THEN
ARG0(1)=2.D00*PI+DIMAG(LOG(Z1(1)/Z2(1)))
ELSE
ARG0(1)=DIMAG(LOG(Z1(1)/Z2(1)))
END IF
C IF(DIMAG(Z1(2)/Z2(2)).LT.0.d00) THEN
C ARG0(2)=2.D00*PI+DATAN2(DIMAG(Z1(2)/Z2(2)),DREAL(Z1(2)/Z2(2)))
C ELSE
C ARG0(2)=DATAN2(DIMAG(Z1(2)/Z2(2)),DREAL(Z1(2)/Z2(2)))
C END IF
IF(DIMAG(LOG(Z1(2)/Z2(2))).LT.0.d00) THEN
ARG0(2)=2.D00*PI+DIMAG(LOG(Z1(2)/Z2(2)))
ELSE
ARG0(2)=DIMAG(LOG(Z1(2)/Z2(2)))
END IF
C
DO I=1,2
DO J=1,2
C(I,J)=0.D00
END DO
END DO
DO I=1,2
DO J=1,2
DO L=1,2
C(I,J)=C(I,J)+DREAL(QP(J,I,L)*(DLOG(CMOD(Z1(L))/CMOD(Z2(L)))
@ +IC*ARG0(L)))
END DO
END DO
END DO
RETURN
END
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
