// JOB
// FOR
*LIST ALL
*IOCS(2501 READER, 1403 PRINTER)
*IOCS(DISK)
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
C LA DESCRIPCION DE UNO DE LOS PROGRAMAS QUE SE EMPLEARON PARA
C CALCULAR EL ESPECTRO DE FRECUENCIAS Y LOS MODOS NORMALES DE
C VIBRACION ES LLAMADO PENTA
C PENTA ES UN PROGRAMA PRINCIPAL QUE NOS CALCULA LOS EIGENVALORES
C Y EIGENVECTORES DE UNA MATRIZ PENTADIAGONAL A LA CUAL SE LE
C PUEDEN HACER MODIFICACIONES EN LAS ESQUINAS DE SUS DIAGONALES
C INSTRUCCIONES DE ENTRADA QUE SIRVEN PARA DEFINIR ARREGLOS,
C NOMBRES DE VARIABLES Y PARA DEFINIR ARREGLOS EN EL DISCO
C
DIMENSION Z(32)
COMMON N,LI,LO,XO,XN
COMMON D(30),C(30),B(30),A(30),E(30)
COMMON Q(4,4),V(31)
DEFINE FILE 10(57,120,U,Kl)
DEFINE FILE 101(41,96,U,K01)
DEFINE FILE 102(41,96,U,K02)
DEFINE FILE 103(41,96,U,K03)
DEFINE FILE 104(41,96,U,K04)
DEFINE FILE 105(41,96,U,K05)
DEFINE FILE 106(41,96,U,K06)
DEFINE FILE 107(41,96,U,K07)
DEFINE FILE 108(41,96,U,K08)
DEFINE FILE 109(41,96,U,K09)
DEFINE FILE 110(41,96,U,K10)
DEFINE FILE 111(41,96,U,Kll)
DEFINE FILE 112(41,96,U,K12)
DEFINE FILE 113(41,96,U,K13)
DEFINE FILE 114(41,96,U,K14)
DEFINE FILE 115(41,96,U,K15)
DEFINE FILE 116(41,96,U,K16)
DEFINE FILE 117(41,96,U,K17)
DEFINE FILE 118(41,96,U,K18)
DEFINE FILE 119(41,96,U,K19)
DEFINE FILE 120(41,96,U,K20)
C
C INSTRUCCIONES DE ENTRADA PARA LEER LA DIMENSION (N) DE LA MATRIZ
C QUE DESEAMOS ANALIZAR
C
LI = 8
LO = 5
1 READ (LI,20) N
20 FORMAT (I2)
IF (N) 2,2,3
2 CALL EXIT
3 WRITE (L0,300) N
300 FORMAT ('IN='I2)
C
C GENERACION DE LOS ELEMENTOS MATRICIALES.- LLENAMOS LA MATRIZ
C TENIENDO EN CUENTA SU FORMA PENTADIAGONAL (DIAGONALES A,B,C,D,E)
C LO CUAL NOS PERMITE UTILIZAR MENOS MEMORIA DE LA MAQUINA, ESTA
C PARTE ES LA QUE GENERALMENTE DIFIERE DE MODELO A MODELO
C
X = -0.8
X = 0.0
X = 0.05
DO 30 K=1,41
X1 = X+1.0
X2 = 1.0/X1
X3 = 1.0-X2
DO 10 J=1,N
A(J) = X2
B(J) = X3
C(J) = -2*(X2+X3)
D(J) = X3
10 E(J) = X2
C
C PENII CALCULA LOS LIMITES DE GERSHORIN DE LOS EIGENVALORES Y
C FORMA LOS COEFICIENTES USADOS EN LA RELACION DE RECURSION COMO
C ELEMENTOS DE LAS MATRICES DE TRANSFERENCIA (T)
C
CALL PENII
C
C PENRO LOCALIZA LAS RAICES DE LA ECUACION CARACTERISTICA,
C COLOCANDO A ELLAS EN UN CIERTO ARREGLO (V)
C
CALL PENRO
C
C PENGR IMPRIME ESTRELLAS QUE REPRESENTAN LOS VALORES NUMERICOS
C DE LOS EIGENVALORES
C
CALL PENGR (V,8.0)
C
C CALCULO DE LOS EIGENVECTORES
C
DO 12 I=1,N
C
C PENEV CALCULA LAS COMPONENTES DEL I-ESIMO EIGENVECTOR, DEJANDO
C LAS COMPONENTES EN EL ARREGLO Z
C
CALL PENEV (Z,I)
C
C PENNO NORMALIZA LAS PRIMERAS N COMPONENTES DEL VECTOR DE
C ARGUMENTO Z
CALL PENNO (Z,N)
C
C PENST(Z,N,M,K,L) ALMACENA Y RECUPERA LOS VECTORES DE EL DISCO
C DE ACUERDO CON LA OPCION L, EL ARGUMENTO ES UN VECTOR Z DE
C DIMENSIONES N,K ES EL ARCHIVO DEL DISCO QUE SERA EMPLEADO Y M
C ES EL NUMERO DE REGISTRO
C
12 CALL PENST (Z,N,100+I,K,3)
C
C EN ESTA PARTE GRAFICAMOS LOS MODOS NORMALES DE VIBRACION LA
C OPCION ES QUE PUEDEN SER GRAFICADOS CON ESTRELLAS Y UNIDOS
C MANUALMENTE O GUARDADOS EN ALGUN REGISTRO PARA QUE DESPUES
C SEAN UNIDOS AUTOMATICAMENTE
C
30 X = X + 0.05
DO 32 I=1,N
C
C PENPG(Z,N,M,K,L) HACE UNA PAGINA DE GRAFICAS DE LOS EIGENVEC-
c TORES, Z ES EL VECTOR QUE CONTIENE LAS COMPONENTES (NORMALIZADAS)
C DE EL EIGENVECTOR DE DIMENSION N, M ES LO QUE DESEAMOS SACAR, K
C ES EL NUMERO DE ARCHIVO EN QUE LOS PUNTOS VAN A SER INTRODUCIDOS
C Y L ES LA OPCION
C
CALL PEMPG (Z,N,0,10,1)
CALL PENPG (Z,N,0,10,2)
CALL PENST (Z,N,100+I,41,4)
DO 31 K=1,41
CALL PENST (Z,N,100+I,K,5)
31 CALL PENPG (Z,N,K,10,3)
32 CALL PENPG (V,N,I,10,5)
GO TO 1
END
//EJECT
// JOB
// FOR
*LIST ALL
*IOCS(2501 READER, 1403 PRINTER)
*IOCS(DISK)
*ONE WORD INTEGERS
*EXTENDED PRECISION
*NAME JOC
C JOC ES EL PROGRAMA PRINCIPAL DE UN CONJUNTO DE PROGRAMAS, DESIGNADOS
C PARA EVALUAR LOS EIGENVALORES Y EIGENVECTORES DE UNA MATRIZ PENTADIAGONAL
C
COMMON N,LI,LO,XO,XN
COMMON P(30,5)
COMMON Q(4,4),V(31)
WRITE (3,300)
300 FORMAT('1'20X' DIAGONALIZACION DE UNA MATRIZ PENTADIAGONAL '/
1 30X' 14 DE FEBRERO DE 1970 '/)
1 CALL PENRM
CALL PENII
CALL PENRO
CALL PENWM
WRITE (3,350)
350 FORMAT('1')
GO TO 1
END
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENCP
C PENCP IMPRIME EL VALOR DEL POLINOMIO CARACTERISTICO EN 100
C PUNTOS ASI COMO LOS ELEMENTOS DEL MENOR QUE LO DEFINE Y LA FASE
C DE SUS VECTORES COLUMNA
C
COMMON N,LI,LO,XO,XN
COMMON P(30,5)
COMMON Q(4,4),V(31)
EQUIVALENTE (Q11,Q(1,1)),(Q12,Q(1,2)),(Q21,Q(2,1)),(Q22,Q(2,2))
WRITE (3,300)
300 FORMAT('1 POLINOMIO CARACTERISTICO DE UNA MATRIZ PENTADIAOONAL'/)
X = XO - 0.5
DX = (XN-XO+1.0)*100.0
DO 2 I=1,100
CALL PENQ(X,Y,N)
P1 = ATAN(Q12/Q11)
P2 = ATAN(Q22/Q21)
WRITE (3,310) Q,Y,Q11,Q12,PI,Q21,Q22,P2
310 FORMAT ('0X='F8.3,3X,'Q(X)='E16.9,
* 5X'Q11='F15.6', Q12='F16.6,10X,'PHI='F12.8/
* 40X'Q21='F15.6', Q22='F15.6,10X,'PHI='F12.8)
2 X = X+DX
WRITE (3,350)
350 FORMAT ('1')
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENEV (Z,K)
C
C PENEV(Z.K) CALCULA LAS COMPONENTES DEL K-ESIMO EIGENVECTOR, DEJANDO
C LOS ELEMENTOS EN EL ARREGLO Z
DIMENSION S(32),T(34),Z(32)
COMMON N,LI,LO,XO,XN
COMMON A(30),B(30),C(30),D(30),B(30)
COMMON Q(4,4),V(31)
EQUIVALENCE (S(1),T(3))
EQUIVALENCE (Q11,Q(1,1)),(Q12,Q(1,2))
CALL PENQK (V(K),Y,N)
T(1) = 0.0
T(2) = 0.0
T(3) = Qll
T(4) = -Q12
DO 10 J=1,N
10 S(J+2)=A(J)*T(J+3)+(B(J)+V(K)*E(J))*T(J+2)+C(J)*T(J+1)+D(J)*T(J)
N2 = N+2
DO 30 I=1,N2
30 Z(I) = S(I)
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENGE (V)
C
C PENGE(V) SE EMPLEA PARA SACAR LOS EIGENVECTORES DEL DISCO LOS UNIFICA
C Y GRAFICA, POR LO REGULAR ES UNA PARTE DEL PROGRAMA PRINCIPAL, ESA ES
C LA RAZON POR LA CUAL NO APARECE EN ALGUNOS DE ELLOS
C
DIMENSION V(1),Z(31)
COMMON N
DO 32 I=1,N
CALL PENPG (Z,N,0,10.1)
CALL PENPG (Z,N,0,10,2)
CALL PENST (Z,N,100+I,41,4)
DO 31 K=1,41
CALL PENST (Z,N,100+I,K,5)
31 CALL PENPG (Z,N,K,10,3)
32 CALL PENPG (V,N,I,10,5)
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENGM (A1,A2,L)
C
C PENGM (A1,A2,L) SIMULA LA ENTRADA DE TARJETAS DE DATOS PARA UNA
C MATRIZ PENTADIAGONAL CON LA PRIMERA DIAGONAL EXTERNA IGUAL A A1 Y LA
C SEGUNDA EXTERNA IGUAL A A2, EXCEPTO PARA LAS OPCIONES L=1,2,3,4 LAS
C TERMINACIONES SON
C L=1 AMBOS EXTREMOS FIJOS
C L=2 EXTREMO IZQUIERDO FIJO, DERECHO LIBRE
C L=3 EXTREMO IZQUIERDO LIBRE, DERECHO FIJO
C L=4 AMBOS EXTREMOS LIBRES
C
DIMENSION A(30),B(30),C(30),D(30),E(30)
COMMON N,LI,LO,XO,XN
COMMON P(30,5)
EQUIVALENCE (A(1),P(1,4)),(B(1),P(1,3)),(C(1),P(1,2))
EQUIVALENCE (D(1),P(1,1)),(E(1),P(1,5))
DO 1 I=1,N
A(I) = A2
B(I) = Al
C(I) = -2.0*(A1+A2)
D(I) = Al
E(I) = A2
1 CONTINUE
GOTO (10,20,30,40) ,L
10 RETURN
20 C(N-l) = C(N-1)+A2
C(N) = C(N)+A1+A2
RETURN
30 C(1) = C(1)+A2+A1
C(2) = C(2) +A2
RETURN
40 C(1) = C(1)+A1+A2
C(N) = C(1)
C(2) = C(2)+A2
C(N-1)=C(2)
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENGR (Z,W)
C
C PENGR IMPRIME ESTRELLAS QUE REPRESENTAN LOS VALORES NUMERICOS DE LOS
C EIGENVALORES Y ADEMAS COLOCA UNA '@' EN EL PUNTO DONDE SE INTERSECTAN
C DOS CURVAS
C
DIMENSION Z(1),III(120)
COMMON N,LI,LO
S=100.0/W
DO 11 I=1,120
11 III (I)=16448
DO 12 I=1,101,10
12 III(I)=-14016
DO 27 I=1,N
L=IFIX(-Z(I)*S)+1
IF(L) 27,27,21
21 IF (L-120) 22,22,27
22 IF (III(L)-23616) 23,24,23
23 III(L) = 23616
GO TO 27
24 III(L)=20544
27 CONTINUE
WRITE (L0.300) III
300 FORMAT (1X,120A1)
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENII
C
C PENNI CALCULA LOS LIMITES DE GERSHGORIN DE LOS EIGENVALORES Y
C FORMA LOS COEFICIENTES USADOS EN LAS RELACIONES DE RECURSION COMO
C ELEMENTOS DE LAS Q'S
C
DIMENSION A(30),B(30),C(30),D(30),E(30)
COMMON N,LI,LO,XO,XN
COMMON P(30,5)
EQUIVALNECE (A(1),P(1,4)),(B(1),P(1,3) ),(C(1),P(1,2))
EQUIVALENCE (D(1),P(1,1)),(E(1),P(1,5))
C
C ESTABLECE LOS LIMITES DE GERSHGORIN PARA LOS EIGENVALORES
C
R=ABS(D(1))+ABS(E(1))
XO=C(1)-R
XN=C(1)+B
R =ABS(B(2))+ABS(D(2))+ABS(E(2))
XO=PENMI(C(2)-R,XO)
XN=PENMA(C(2)+R,XN)
N2=N-2
DO 4 I=3,N2
R =ABS(A(I))+ABS(B(I))+ABS(D(I))+ABS(E(I))
XO=PENMI(C(I)-R,XO)
4 XN=PENMA(C(I)+R,XN)
R =ABS(A(N-1))+ABS(B(N-1))+ABS(D(N-1))
XO=PENMI(C(N-1)-R,XO)
XN=PENMA(C(N-1)+R,XN)
R =ABS(A(N))+ABS(B(N))
XO=PENMI(C(N)-R,XO)
XN=PENMA(C(N)+R,XN)
C
C FORMAMOS LOS COEFICIENTES QUE SERAN USADOS EN LAS RELACIONES DE
C RECURSION
C
DO 20 I=1,N
R = 1.0/E(I)
DO 16 J=1,4
16 P(I,J)=-R*P(I,J)
20 E(I)=R
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
FUNCTION PENMA(X.Y)
C
C PENMA(X,Y) SELECCIONA EL MAXIMO DE X Y Y
C
IF (X-Y) 1,2,2
1 PENMA=Y
RETURN
2 PENMA=X
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
FUNCTION PENMI(X,Y)
C
C PENMI(X,Y) SELECCIONA EL MINIMO DE X Y Y
C
IF (X-Y) 1,1,2
1 PENMI=X
RETURN
2 PENMI=Y
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENNO (Z,N)
C
C PENNO(Z,N) NORMALIZA LAS PRIMERAS COMPONENTES DE EL VECTOR DE
C ARGUMENTO Z
C
DIMENSION Z(1)
A = 0.0
DO 5 I=1,N
5 A = A+Z(I)**2
A = 1.0/SQRT(A)
DO 7 I=1,N
7 Z(I)=A*Z(I)
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENPG (Z,N,M,K,L)
C
C PENPG(Z,N,M,K,L) HACE UNA PAGINA DE GRAFICAS DE LOS EIGENVECTORES,
C Z ES EL VECTOR QUE CONTIENE LAS COMPONENTES (NORMALIZADAS) DEL
C EIGENVECTOR DE DIMENSION N, M ES EL QUE DESEAMOS SACAR, K ES EL
C NUMERO DE ARCHIVO EN EL QUE LOS PUNTOS VAN A SER INTRODUCIDOS Y
C L ES LA OPCION DESEADA
C L=1 GENERA UNA PAGINA EN BLANCO
C L=2 INVIERTE LAS COMPONENTES DEL EIGENVECTOR
C L=3 IMPRIME LA PAGINA TERMINADA
C
DIMENSION Z(1),JJJ(120)
COMMON NN,LI,LO
GOTO (10,20,30,40,50) ,L
10 CONTINUE
DO 11 I=1,120
11 JJJ(I)=16448
DO 12 I=1,57
12 WRITE (K'I) JJJ
RETURN
20 CONTINUE
DO 22 J=1,41
READ (K'50-J) JJJ
DO 21 I=1,N
JJ=7*I+J
21 JJJ(JJ)=19264
22 WRITE (K'50-J) JJJ
RETURN
30 CONTINUE
DO 31 I=1,N
31 CALL PENPL (7*I+M,50 -(IFIX(10.0*Z(I))+M),K,23616)
RETURN
40 CONTINUE
WRITE(K'M) (Z(J),J=1,N)
RETURN
C
C CUANDO L=5 LOS PARAMETROS SON USADOS COMO SIGUE:
C Z ES EL VECTOR EIGENVALOR
C M LA COMPONENTE DEL VECTOR FRECUENCIA QUE SERA USADO
C K EL ARCHIVO QUE SERA IMPRIMIDO
C
353 FORMAT(' MODO NORMAL NUMERO',I5.5X,'EIGENVALOR',F15.6/)
50 CONTINUE
WRITE (LO,352)
WRITE (LO,353) M,Z(M)
WRITE (LO,350)
DO 51 I=1,57
READ (K'I) JJJ
WRITE (LO,351) JJJ
51 CONTINUE
WRITE (LO,350)
RETURN
350 FORMAT (1X,12('P E N T A '))
351 FORMAT (1X,120A1)
352 FORMAT (1H1)
353 FORMAT(' MODO NORMAL NUMERO',I5,5X,'EIGENVALOR',F15.6/)
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENPL (IA,IB,K,ICH)
C
C PENPL(IA,IB,K,ICH) INTRODUCE EL CARACTER ICH EN EL ARCHIVO DEL
C DISCO K CON LA COLUMNA IA DEL NUMERO DE REGISTRO IB
C
DIMENSION III(120)
IF (IA) 1,1,2
1 IA=1
2 IF (120-IA) 3,3,4
3 IA=120
4 IF (IB) 5,5,6
5 IB=1
6 IF (57-IB) 7,8,8
7 IB=57
8 READ (K'IB) III
III(IA)=ICH
WRITE (K'IB) III
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENQK (X,Y,K)
C
C PENQK(X,Y,K) FORMA EL PRODUCTO DE LAS PRIMERAS K DE LAS Q-MATRICES
C EVALUANDO EL PRODUCTO RECURSIVAMENTE CON EL EIGENVALOR X, Y DEPOSITA
C EL VALOR DEL MENOR SUPERIOR IZQUIERDO 2X2 DE EL PRODUCTO EN Y. SOLAMEN-
C TE LAS PRIMERAS DOS COLUMNAS DE Q SON EVALUADAS ASI COMO SOLAMENTE
C ELLAS ENTRAN DENTRO DE LA EVALUACION DEL POLINOMIO CARACTERISTICO
C
COMMON N,LI,LO,XO,XN
COMMON A(30),B(30),C(30),D(30),E(30)
COMMON Q(4,4),V(31)
EQUIVALENCE (Qll,Q(1,1)),(Q12,Q(1,2)),(Q21,Q(2,1)),(Q22,Q(2,2))
EQUIVALENCE (Q31,Q(3,1)),(Q32,Q(3,2)),(Q41,Q(4,1)),(Q42,Q(4,2))
IF (K-1) 1,2,3
1 Y = 1.0
RETURN
2 Y=B(1)+X*E(1)
RETURN
3 Q11=A(1)
Q12=B(1)+X*E(1)
Q21=1.0
Q22=0.0
Q31=0.0
Q32=1.0
Q41=0.0
Q42=0.0
W1=Q11
W2=Q12
DO 20 J=2,K
W1=(B(J)+X*E(J))*Q21+C(J)*Q31+D(J)*Q41
W2=(B(J)+X*E(J))*Q22+C(J)*Q32+D(J)*Q42
Q41=Q31
Q42=Q32
Q31=Q21
Q32=Q22
Q21=Q11
Q22=Q12
Q11=W1+A(J)*Q11
Q12=W2+A(J)*Q12
20 CONTINUE
Y=W1*Q22-W2*Q21
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENRM
C
C PENRM LEE LAS TARJETAS DE ACUERDO CON LOS ELEMENTOS PENTADIAGONALES
C DE LA MATRIZ Y FORMA LAS FILAS DE ARRIBA DE LAS Q- MATRICES, N Y
C EPSILON SON ESPECIFICADOS EN LA PRIMERA TARJETA
C
COMMON N,LI,LO,XO,XN
COMMON D(30),C(30),B(30),A(30),E(30)
READ (LI,200) N
200 FORMAT (I2)
IF (N) 1,1,2
1 RETURN
2 WRITE (LO,300) N
300 FORMAT (' N= ',I2/)
WRITE (LO,301)
301 FORMAT('O MATRIZ PENTADIAGONAL QUE SERA DIAGONALIZADA')
DO 3 I=1,N
READ (LI,210) A(I),B(I),C(I),D(I),E(I)
210 FORMAT (5F10.0)
WRITE (LO,310) A(I),B(I),C(I),D(I),E(I)
310 FORMAT (1X,5F15.6)
3 CONTINUE
RETURN
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENRO
C
C PENRO LOCALIZA LAS RAICES DE LA ECUACION CARACTERISTICA, COLOCANDO
C A ELLAS EN EL ARREGLO V
C
DIMENSION J(2),R(2),S(2),T(31),U(32)
COMMON N,LI,LO,XO,XN
COMMON A(30),B(30),C(30),D(30),E(30)
COMMON Q(4,4),V(31)
EQUIVALENCE (Q11,Q(1,1)),(Q12,Q(1,2)),(Q21,Q(2,1),(Q22,Q(2,2))
EQUIVALENCE (J1,J(1)),(J2,J(2)),(R1,R(1)),(R2,R(2))
EQUIVALENCE (S1,S(1)),(S2,S(2))
SIZ(X) = ABS(X)-EPS
WID(X,Y) = ABS(Y- X)- DEL
QM(X1,X2,X3,Y1,Y2,Y3)=(0.5*(Y1*(X2-X3)*(X2+X3)+Y2*(X3-X1)*(X3+X1)
* +Y3*(X1-X2)*(X1+X2)))/(Y1*(X2-X3)+Y2*(X3-X1)
* +Y3*(X1-X2))
RS(X,Y) = SIGN(l.0,X)-SIGN(1.0,Y)
PHI = 0.15
TAU = l.OE-2
DEL = (1.0E-5*ABS(XN-XO))/N
U(1) = XO-DEL
C
C PRODUCIMOS LA PRIMERA RAIZ CUYO VALOR ES EL ELEMENTO A11 DE LA MATRIZ
C
V(1) = -B(1)/E(1)
C
C LOS INTERVALOS CERRADOS QUE CONTIENEN LAS RAICES SON DETERMINADOS
C RECURSIVAMENTE
C
DO 500 K=2,N
C
C LOS CEROS DE LOS POLIMOMIOS ANTERIORES DETERMINAN INTERVALOS PARA
C ESTE CICLO, ASI COMO EL ESTABLECIMIENTO DEL FACTOR DE ESCALA QUE
C SERA EMPLEADO
C
CALL PENQK (U(1),T(1),K)
DO 1 I=2,K
U(I) = V(I-1)
1 CALL PENQK (U(I),T(I),K)
U(K+1) = XN
CALL PENQK (U(K+1),T(K+1),K)
U(N+2) = XN+DEL
C
C UNA RAIZ TIENE QUE SERA LOCALIZADA PARA CADA K INTERVALO
C
DO 400 I=1,K
X1 = U(I)
Y1 = T(I)
X2 = U(I+1)
Y2 = T(I+1)
X3 = U(I+2)
IF (RS(Y1,Y2)) 50,5,50
C
C SI EL INTERVALO ES EXTREMADAMEMTE PEQUENO, SUPONDREMOS QUE LA RAIZ
C QUE ESTA CONTENIDA SE ENCUENTRA EN SU PUNTO MEDIO
C
5 IF (WID(X1,X2)) 6,6,7
6 V(I) = 0.5*(X1+X2)
GO TO 400
7 R1 = X2-PHI*(X2-X1)
R2 = X2+PHI*(X3-X2)
DO 25 M=1,3
DO 15 L=1,2
CALL PENQK (R(L),S(L),K)
IF (RS(Y1, S(L))) 12,15,12
12 U(I+1) = R(L)
T(I+1) = S(L)
X2 = R(L)
Y2 = S(L)
GO TO 40
15 CONTINUE
X = QM(Rl,X2,R2,Sl,Y2,S2)
CALL PENQK (X,Y,K)
IF (Xl-X) 16,16,23
16 IF (X-X3) 17,17,23
17 IF (RS(Y1,Y)) 19,20,19
19 U(I+1) = X
T(I+1) = Y
X2 = X
Y2 = Y
GO TO 40
20 X2 = X
Y2 = Y
22 R1 = 0.5*(R1+X)
R2 = 0.5*(X+R2)
GO TO 25
23 X = X2
GO TO 22
25 CONTINUE
30 V(I) = X-TAU*(X-X1)
V(I+1)= X+TAU*(X3-X)
I = I+1
GO TO 400
40 IF (U(I+1)-U(I+2)) 50,50,41
41 U(I+2) = U(I+1)
T(I+2) = T(I+1)
GO TO 50
C
C ***************************************
C * RUTINA QUE LOCALIZA RAICES *
C ***************************************
C
50 A1 = (X1+X1+X2)/3.0
A2 = (X1+X2+X2)/3.0
CALL PENQK (A1,Z1,K)
CALL PENQK (A2,Z2,K)
EPS = 0.5E-6*(ABS(Z1)+ABS(Z2))
C
C EMPLEANDO EL METODO DE LA "REGLA FALSA" SE ENCUENTRA NORMALMENTE QUE
C EL VALOR DE LA FUNCION EN UN PUNTO EXTREMO ES MUCHO MENOR QUE EN EL
C OTRO CUAMDO ES PEQUENO EL PUNTO EXTREMO PUEDE SER TOMADO COMO UNA
C RAIZ CUANDO LA RAZON DE LOS VALORES EXTREMOS ES MENOR QUE LA PRECISION
C DE LA COMPUTADORA, LA NUEVA APROXIMACION SERA NO CONFIABLE Y EL
C PUNTO EXTREMO MAYOR TIENE QUE SER MOVIDO.
DO 200 M=1,30
X = 0.5*(X1+X2)
CALL PENQK (X,Y,K)
IF (RS(Y,Y1)) 106,105,106
105 X1 = X
Y1 = Y
GO TO 107
106 X2 = X
Y2 = Y
107 IF (WID(X1,X2)) 108,108,110
108 V(I) = 0.5*(X1+X2)
GO TO 400
110 IF (ABS(Y1/Y2+Y2/Y1)-2.0E4) 125,200,200
125 Sl = X1
S2 = X2
C
C DOS CICLOS DE LA REGLA FALSA SE LLEVAN A CABO, AUNQUE NO ES LOGICA-
C MENTE NECESARIO VERIFICAR QUE LA NUEVA ESTIMACION SE ENCUENTRE DENTRO
C DEL INTERVALO DE RAIZ, LA PERDIDA DE PRECISION PUEDE ALGUNAS VECES
C ARROJAR LA INTER-SECCION FUERA.
C
DO 135 L=1,2
R(L) = (X1*Y2-X2*Y1)/(Y2-Yl)
IF (Xl-R(L)) 129,200,200
129 IF (R(L)-X2) 130,200,200
130 CALL PENQK (R(L),Y,K)
IF (SIZ(Y)) 131,131,132
131 V(I)=R(L)
GO TO 400
132 IF (RS(Y,Y1)) 134,133,134
133 XI = R(L)
Yl = Y
J(L) = -1
GO TO 135
134 X2 = R(L)
Y2 = Y
J(L) = 1
135 CONTINUE
C
C EL PROCESO DE AITKEN DEL CUADRADO DE LA DELTA PROPORCIONA LAS ULTIMAS
C DOS APROXIMACIONES TRASLADADAS EN EL MISMO PUNTO EXTREMO Y PROPORCIONA
C LA NUEVA ESTIMACION LOCALIZADA DENTRO DEL INTERVALO
C
IF (J1+J2) 141,200,142
141 RO = S1
GO TO 150
142 RO = S2
150 X = (RO*R2-R1*R1)/(RO-R1-R1+R2)
IF (Xl-X) 151,200,200
151 IF (X-X2) 152,200,200
152 CALL PEMQK (X,Y,K)
IF (SIZ(Y)) 153,153,154
153 V(I) = X
GO TO 400
154 IF (RS(Y,Y1)) 156,155,156
155 X1 = X
Y1 = Y
GO TO 200
156 X2 = X
Y2 = Y
GO TO 200
200 CONTINUE
400 CONTINUE
500 CONTINUE
RETURN
END
C
// JOB
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENST (Z,N,K,M,L)
C
C PENST(Z,N,K,M,L) ALMACENA Y RECUPERA LOS VECTORES DEL DISCO DE
C ACUERDO CON LA OPCION L, EL ARGUMENTO ES UN VECTOR Z DE DIMENSION
C N, K ES EL ARCHIVO DE EL DISCO QUE SERA EMPLEADO Y M ES EL NUMERO
C DE REGISTRO
C L=1
C L=2
C L=3 ALMACENA EL VECTOR EN EL DISCO
C L=4 REVISA TODOS LOS VECTORES PARA LA MAYOR CONTINUIDAD POSIBLE
C L=5 BUSCA EL VECTOR DEL DISCO
C
DIMENSION X(32),Y(32),Z(32)
GOTO (10,20,30,40,50) ,L
10 CONTINUE
RETURN
20 CONTINUE
RETURN
30 CONTINUE
WRITE (K'M) (Z(I),I=1,N)
RETURN
40 CONTINUE
READ (K'l) (X(I),I=1,N)
DO 46 J=2,M
READ (K'J) (Y(I),I=1,N)
A = 0.0
B = 0.0
DO 41 I=1,N
A = A+(Y(I)-X(I))**2/((I+1)*(N-I+1) )
B = B+(Y(I)+X(I))**2/((I+1)*(N-I+1))
41 CONTINUE
IF (A-B) 44,44,42
42 DO 43 I=1,N
43 Y(I) = -Y(I)
44 DO 45 I=1,N
45 X(I) = Y(I)
46 WRITE (K'J) (Y(I),I=1,N)
RETURN
50 CONTINUE
READ (K'M) (Z(I),I=1,N)
RETURN
END
C
// JOB
// FOR
*LIST ALL
*IOCS(2501 READER, 1403 PRINTER)
*IOCS(DISK)
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENW1 (XI,Y,L)
C
C PENW1(XI,Y,RHO,L) GRAFICA LOS PUNTOS EN EL PLANO XI-Y COMO FUNCIONES
C DE LOS DOS COSENOS DE LOS NUMEROS DE ONDA, VARIAS OPCIONES EXISTEN DE
C ACUERDO CON LA OPCION L
C L=1 LIMPIA LA IMAGEN DE LA PAGINA E INSERTA LA FRONTERA C=1
C L=2 INSERTA EL PUNTO DE ACUERDO CON LA REL. DE DISPERSION MONOATOMICA
C L=3 INSERTA EL PUNTO DE ACUERDO CON LA REL. DE DISPERSION DIATOMICA
C L=4 IMPRIME LA PAGINA COMPLETA
C
DIMENSION III(120)
COMMON N,LI,LO,XO,XN
IX(X) =92+IFIX(26.0*X)
IY(Y) =18-IFIX(15.0*Y)
GO TO (10,20.40),L
10 DO 12 I=1,120
12 III(I) = 16448
DO 14 I=1,57
14 WRITE (20'I) III
CALL PENPL (IX( 0.0),IY( 0.0),20,19264)
CALL PENPL (IX( 1.0),IY( 1.0),20,19264)
CALL PENPL (IX( 1.0),IY(-1.0),20,19264)
CALL PENPL (IX(-l.0),IY( 1.0),20,19264)
CALL PENPL (IX(-l.0),IY(-1.0),20,19264)
RETURN
C
C LOS ARGUMENTOS PARA PENGW DEL MODELO MONOATOMICO SON
C XI = -LAMBDA/(4.0*(A1+A2))
C Y = A1/(4.0*A2)
C RHO=1.0
C
20 B = Y+1.0
AC = XI*(1.0+4.0*Y)
IF (B*B-AC) 22,23,24
22 CONTINUE
SQ = SQRT(AC-B*B)
C1 = -Y+SQ
C2 = -Y-SQ
CALL PENPL (IX(C1),IY(C2),20,23616)
RETURN
23 CALL PENPL (IX(B),IY(B),20,23616)
RETURN
24 CONTINUE
SQ = SQRT(B*B-AC)
C1 = -Y+SQ
C2 = -Y-SQ
CALL PENPL (IX(C2),IY(C1),20,23616)
RETURN
40 WRITE (LO,3010)
DO 42 I=1,57
READ (20'I) III
WRITE (LO,3000) III
42 CONTINUE
RETURN
3000 FORMAT (1X,120A1)
3010 FORMAT ('1')
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENW2 (XI, Y, RHO, L)
C
C PENW2(XI,Y,RHO,L) GRAFICA LOS PUNTOS EN EL PLANO XI-Y COMO
C FUNCIONES DE LOS DOS COSENOS DE LOS NUMEROS DE ONDA, VARIAS OPCIONES
C EXISTEN DE ACUERDO CON LA OPCION L
C L=1 LIMPIA LA IMAGEN DE LA PAGINA E INSERTA LA FRONTERA C=1
C L=2 INSERTA EL PUNTO DE ACUERDO CON LA REL. DE DISPERSION MONOATOMICA
C L=3 INSERTA EL PUNTO DE ACUERDO CON LA REL. DE DISPERSION DIATOMICA
C L=4 IMPRIME LA PAGINA COMPLETA
DIMENSION III(120)
COMMON N,LI,LO,XO,XN
IX(X) = 92+IFIX(26.0*X)
IY(Y) = 40-IFIX(15.0*X)
GO TO (10,30,40),L
10 DO 12 I=1,120
12 III(I) = 16448
DO 14 I=1,57
14 WRITE (20'I) III
CALL PENPL (IX( 0.0),IY( 0.0),20,19264)
CALL PENPL (IX( 1.0),IY( 1.0),20,19264)
CALL PENPL (IX( 1.0),IY(-1.0),20,19264)
CALL PENPL (IX(-l.0),IY( 1.0),20,19264)
CALL PENPL (IX(-l.0),IY(-1.0),20,19264)
RETURN
C
C LOS ARGUMENTOS PARA PENGW DEL MODELO DIATOMICA SON
C XI = LAMBDA/A2
C Y = A1/A2
C RHO= M PEQUENA/ M GRANDE
C
30 CRHO = SQRT(RHO)
CRHO = CRHO+1.0/CRHO
B = CRHO*XI+(Y+2.0)**2
AC = 4.0*(XI*XI+2.0*CRHO*XI*(Y+1)+2.O*((Y+2.0)**2-2.0))
IF (B*B-AC) 32,33,34
32 SQ = SQRT(AC-B*B)
C1 = 0.25*(B+SQ)
C2 = 0.25*(B-SQ)
CALL PENPL (IX(C1),IY((C2),20,23616)
RETURN
33 CALL PENPL (IX(B),IY(B),20,23616)
RETURN
34 SQ = SQRT(B*B)-AC)
C1 = 0.25*(B+SQ)
C2 = 0.25*(B-SQ)
CALL PENPL (IX(C2),IY(C1),20,23616)
RETURN
40 WRITE (LO,3010)
DO 42 I=1,57
READ (20'I) III
WRITE (LO,3000) III
42 CONTINUE
RETURN
3000 FORMAT (1X,120A1)
3010 FORMAT ('1')
END
C
// EJECT
// FOR
*LIST ALL
*ONE WORD INTEGERS
*EXTENDED PRECISION
C
SUBROUTINE PENWM
C
C PENWM ES UNA SUBRUTINA QUE ESCRIBE LOS VALORES DE LOS EIGENVALORES Y
C SUS EIGENVECTORES ASOCIADOS
C
DIMENSION Z(32)
COMMON N,LI,LO,XO,XN
COMMON P(30,5)
COMMON Q(4,4),V(31)
WRITE (LO,300)
300 FORMAT ('1')
DO 30 I=1,N
CALL PENEV (Z,I)
R1 = Z(N+1)/Z(1)
R2 = Z(N+2)/Z(1)
CALL PENNO (Z,N)
WRITE (LO,310) I,V(I),(Z(J),J=1,N)
310 FORMAT('O EIGENVECTOR NO.',I2,10X,'EIGENVALOR'
1 ,F10.6/(1X,10F12.6))
WRITE (LO,320) R1,R2
310 FORMAT (' CHECA COMPONENTES QOE RESULTARAN CERO ',10X,2E12.4)
30 CONTINUE
RETURN
END