Matlab - Ayuda urgenbte

Aqui esta el codigo pero da un error al ingresar datos en sentido y si alguien me puede ayudar muchas gracias .... [email protected]

%
ntot=input ('\n Numero total de porticos de la estructura :');
fprintf ('\n Codigos para analisis sismico: Sentido X=1 Sentido Y=2');
isismo=input ('\n Ingrese codigo de sentido de analisis sismico :');
rs=zeros(3*NP,3);rb=[0;0;0];
if isismo==1;nx=iejes; ny=ntot-nx;rb(1)=1; var=2;for j=1:NP; rs(j,1)=1; end;
else;ny=iejes; nx=ntot-ny;rb(2)=1;var=1;for j=1:NP; rs(j+NP,2)=1; end; end;
%Submatrices de rigidez: KEE, con inercias gruesas
Kxx=zeros(NP,NP);Kyy=zeros(NP,NP);Kteta=zeros(NP,NP);cero=zeros(NP,NP);
Kxt=zeros(NP,NP);Kyt=zeros(NP,NP);for k=1:NP;identidad(k,k)=1;end;
for i=1:ntot
for k=1:NP
rtet(k,k)=r(i);
end
rteta=rtet*rtet;
ji=NP*(i-1)+1;jf=NP*(i-1)+NP;
if i<=nx
Kxx=Kxx+KLG(ji:jf,1:NP);Kxt=Kxt+KLG(ji:jf,1:NP)*rtet;
Kteta=Kteta+KLG(ji:jf,1:NP)*rteta;
A(ji:ji+NP-1,:)=[identidad cero rtet];
else
Kyy=Kyy+KLG(ji:jf,1:NP);Kyt=Kyt+KLG(ji:jf,1:NP)*rtet;
Kteta=Kteta+KLG(ji:jf,1:NP)*rteta;
A(ji:ji+NP-1,:)=[cero identidad rtet];
end
end
%Matriz de rigidez espacial con 3 grados de libertad por planta
KS=[Kxx cero Kxt;cero Kyy Kyt;Kxt Kyt Kteta];
fprintf ('\n Matriz de rigidez con base empotrada')
KS
%Matriz de masas de superestructura
disty=abs(r(1))+abs(r(nx));distx=abs(r(nx+1))+abs(r(ntot));
for i=1:NP; masaD(i)=pesoD(i)/9.8; masaL(i)=pesoL(i)/9.8;
mas(i)=masaD(i)+masaL(i);end
masa=zeros(NP,NP);for i=1:NP; masa(i,i)=mas(i);end
mj=zeros(NP,NP);for i=1:NP, mj(i,i)=mas(i)*(distx^2+disty^2)/12;end
MS=[masa cero cero;cero masa cero;cero cero mj];
fprintf ('\n Matriz de masa con base empotrada')
MS
% Matriz de amortiguamiento, tipo Wilson y Penzien (1982)
for i=1:3*NP; zeda(i)=seda; end;zeda=zeda';
CS=zeros(3*NP,3*NP);[V,D]=eig(KS,MS);W=sqrt(diag(D));
% Se ordenan las frecuencias y los modos de vibracion de menor a mayor
[Wn,II]=sort(W); for i=1:3*NP; fi(:,i)=V(:,II(i));end;T1s=2*pi/Wn(1);
for i=1:3*NP
ff=fi(:,i);mi=ff'*MS*ff;aux=2*zeda(i)*Wn(i)/mi;
CS=CS+aux.*MS*ff*ff'*MS;
end
fprintf ('\n Matriz de amortiguamiento con base empotrada')
CS
% Matriz de rigidez del sistema de aislamiento en coordenadas de piso
nd=nx*ny;
fprintf ('\nSi los aisladores son circulares e iguales digite 1, caso contrario 2')
cod= input ('\nIndique el codigo, sobre si los aisladores son iguales o diferentes: ');
Kx=0; Ky=0; Kt=0; Kxte=0; Kyte=0;
if cod==1
kais=input ('Rigidez del aislador: ');
for i=1:nd
kx(i)=kais; ky(i)=kais; Kx=Kx+kx(i); Ky=Ky+ky(i);
end
else
for i=1:nd
fprintf ('\nIndique la rigidez para el aislador %i',i)
kx(i)=input ('\nRigidez del aislador: '); Kx=Kx+kx(i);
ky(i)=kx(i); Ky=Ky+ky(i);
end
end
for i=1:nx
Kt=Kt+kx(i)*r(i)*r(i);Kxte=Kxte+kx(i)*r(i);
end
for i=1:ny
Kt=Kt+ky(i)*r(i+nx)*r(i+nx);Kyte=Kyte+ky(i)*r(i+nx);
end
KB=zeros(3); KB(1,1)=Kx; KB(2,2)=Ky; KB(1,3)=Kxte; KB(2,3)=Kyte;KB(3,3)=Kt;
KB(3,1)=KB(1,3); KB(3,2)=KB(2,3);
fprintf ('\n Matriz de rigidez del sistema de aislacion')
KB
% Matriz de masas del sistema de aislamiento en coordenadas de piso
ij=0;M=0;Jm=0;
for i=1:nx
for j=1:ny
ij=ij+1;
d(ij)=sqrt(r(i)*r(i)+r(j+nx)*r(j+nx));
end
end
fprintf ('\nSi todas las masas son iguales, digite 1, caso contrario 2')
icod= input ('\nIndique el codigo, sobre las masas: ');
if icod==1
masa=input ('\nIndique la masa sobre el aislador: ');
for i=1:nd
m(i)=masa;M=M+m(i);Jm=Jm+m(i)*d(i)*d(i);
end
else
for i=1:nd
fprintf ('\nIndique la masa sobre el aislador %i',i)
m(i)=input('\nMasa sobre el aislador: ');
M=M+m(i); Jm=Jm+m(i)*d(i)*d(i);
end
end
MB=zeros(3); MB(1,1)=M;MB(2,2)=M;MB(3,3)=Jm;
fprintf ('\n Matriz de masas en coordenadas de piso')
MB
% Matriz de amortiguamiento del aislamiento en coordenadas de piso
% Tipo Wilson y Penzien (1982)
for i=1:3; zed(i)=sedabase; end;
CB=zeros(3,3);[VV,DD]=eig(KB,MB);WW=sqrt(diag(DD));
% Se ordenan las frecuencias y los modos de vibracion de menor a mayor
[Wwn,JJ]=sort(WW); for i=1:3; fj(:,i)=VV(:,JJ(i));end; T1b=2*pi/Wwn(1);rel=T1s/T1b;
for i=1:3
fff=fj(:,i);mi=fff'*MB*fff; aux=2*zed(i)*Wwn(i)/mi;
CB=CB+aux.*MB*fff*fff'*MB;
end
fprintf ('\n Matriz de amortiguamiento del sistema de aislacion')
CB
fprintf ('\n Tenga paciencia el programa esta en aislamiento de base, se demora unos segundos')
% Respuesta en el tiempo de aisladores de base
Qo=MB*rb*(-1);
[qt,vt,q]=pse3(MB,CB,KB,Qo,sismo,dt,var);
% Correccion de masa
npuntos=length(sismo);MBI=inv(MB);MBBI=inv(MB');
MSC=MS-MS*rs*MBBI*rs'*MS; MSSC=zeros(3*NP,3*NP); for i=1:3*NP; MSSC(i,i)=MSC(i,i); end
for i=1:npuntos-1;
tt(i)=i*dt;
for j=1:3
qpb(j)=vt(i,j);qb(j)=qt(i,j);
end
if i==1
qpb=qpb';qb=qb';
end
ac=MBBI*(CB*qpb+KB*qb);
if isismo==1; acelera(i)=ac(1); else; acelera(i)=ac(2); end
end
fprintf ('\n Tenga paciencia el programa esta en la estructura, se demora unos segundos')
% Respuesta en el tiempo de superestructura.
Qo=MS*rs;
if isismo==1; Qs=Qo(:,1); else; Qs=Qo(:,2); end
[qs,vs,qq]=pse3(MSSC,CS,KS,Qs,acelera,dt,var);
% Respuestas en porticos extremos, en ultimo piso y en primer piso
for i=1:ntot
ji=NP*(i-1)+1;jf=NP*(i-1)+NP;a=A(ji:jf,1:3*NP);
for j=1:npuntos-2;
t(j)=j*dt;
for k=1:3*NP;
dd(k)=qs(j,k);
end
dd=dd';p=a*dd;dd=dd';
if i==1 & isismo==1; p1(j)=p(NP); end
if i==nx & isismo==1; pnx(j)=p(1); end
if i==nx+1 & isismo==2; p2(j)=p(NP); end
if i==ntot & isismo==2; pnt(j)=p(1); end
end
end
hold on
% Dibujo de respuestas maximas en porticos exteriores en ultimo y 1 piso
%if isismo==1
% p1=p1';pnx=pnx';
% subplot (2,1,1);plot(t,p1); ylabel ('Desplazamiento (m)');title('Desplazamiento ultimo piso. Portico 1');
% subplot (2,1,2); plot(t,pnx); xlabel ('Tiempo (s)');ylabel ('Desplazamiento (m)');
%title ('Desplazamiento primer piso. Portico Extremo en sentido X');
%else
% p2=p2';pnt=pnt';
% subplot (2,1,1);plot (t,p2);ylabel ('Desplazamiento (m)'); title('Desplazamiento ultimo piso. Portico A');
% subplot (2,1,2);plot (t,pnt);xlabel ('Tiempo (s)');ylabel ('Desplazamiento (m)');
% title ('Desplazamiento primer piso. Portico Extremo en sentido Y');
%end
% Dibujo de respuestas maximas de la base y de ultimo piso
plot (tt,q);
plot (t,p1,'r');
%---fin


********************************************

y el error q me despliega es :

Tenga paciencia el programa esta en la estructura, se demora unos segundos??? Undefined function or variable 'p1'.

Error in ==> C:\MATLAB6p5\work\masacorregidaaislamiento.m
On line 219 ==> plot (t,p1,'r');

>>

Podrias decirnos que valor toma NP ?,