I have to extract the values of "p" variable which is calculated in 'model' function. can anyone help?
clear;clc;ifp = fopen('loma.txt','r'); % capture and discard the first two lines of the file
for i=1:46 R = fgetl(ifp); % delete redundant lines
end % load the acceleration data 90 degree
data90 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);ifp = fopen('loma.txt','r'); % capture and discard the first two lines of the filefor i=1:1646 F = fgetl(ifp); % delete redundant linesend % load the acceleration data 0 degree
data0 = fscanf(ifp,'%f')/100; % m/s/s% close the input filefclose(ifp);dt = 0.02;npoint=length(data90);time=0:dt:(npoint-1)*dt;t=time';kx=11*10^4; %N/m
ky=8*10^4; %N/mfyxp=550*10^0; %N
fyxn=650*10^0; %Nfyyp=440*10^0; %Nfyyn=520*10^0; %Nalphax=0.01;alphay=0.01;zux0p=fyxp/kx;zuy0p=fyyp/ky;zux0n=fyxn/kx;zuy0n=fyyn/ky;Z0(1)=0;Z0(2)=0;Z0(3)=0;Z0(4)=0;Z0(5)=0;Z0(6)=0;[t,Z]=ode23s(@model,t,Z0,[],data0,data90,t);qx=alphax*kx*Z(:,1)+(1-alphax)*ky*Z(:,5);qy=alphay*ky*Z(:,3)+(1-alphay)*zuy0p/zux0p*ky*Z(:,6);function [ydot]=model(t,Y,ax2,ay2,t1)ax = interp1(t1,ax2,t);ay = interp1(t1,ay2,t);kx=11*10^4; %N/mky=8*10^4; %N/mfyxp=550*10^0; %Nfyxn=650*10^0; %Nfyyp=440*10^0; %Nfyyn=520*10^0; %Nalphax=0.01;alphay=0.01;beta=0.7;dt = 0.02;eta=1;nu=1;m=50;xi=5/100;s=8;sigma=.2;mu=0.1;zux0p=fyxp/kx;zuy0p=fyyp/ky;zux0n=fyxn/kx;zuy0n=fyyn/ky;zu=((1+sign(Y(5)))*abs(zux0p)+(1-sign(Y(5)))*abs(zux0n))/2/nu;Cx=Y(5)/zu*(1+beta*(sign(Y(2)*Y(5))-1));Cy=Y(6)/zu*(1+beta*(sign(Y(4)*Y(6))-1));px=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(2))*Y(5)/zu-mu)/sigma)^2);py=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(4))*Y(6)/zu-mu)/sigma)^2);H11=eta+(1-Y(5)/zu*Cx)*px;H12=-Y(5)/zu*Cy*py;H21=-Y(6)/zu*Cx*px;H22=eta+(1-Y(6)/zu*Cy)*py;G1=(1-Y(5)/zu*Cx)*Y(2)-Y(5)*zux0p/zu/zuy0p*Cy*Y(4);G2=-Y(6)/zu*Cx*Y(2)+(1-Y(6)/zu*Cy)*zux0p/zuy0p*Y(4);Ydot(5)=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);zdashx=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);Ydot(6)=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);zdashy=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);Ydot(1)=Y(2);Ydot(3)=Y(4);qx=alphax*kx*Y(1)+(1-alphax)*ky*Y(5);qy=alphay*ky*Y(3)+(1-alphay)*zuy0p/zux0p*ky*Y(6);Ydot(2)=-ax-2*xi*sqrt(kx/m)*Y(2)-qx/m;Ydot(4)=-ay-2*xi*sqrt(ky/m)*Y(2)-qy/m;p=1/zu^2*(Y(5)*(Y(2)-mu*zdashx)+Y(6).*(zux0p/zuy0p*Y(4)-mu*zdashy))% v=cumtrapz(dt,p);
% Energy=(1-alphax)/2.*v;
ydot=[Ydot(1);Ydot(2);Ydot(3);Ydot(4);Ydot(5);Ydot(6)];
Best Answer