Hello, I want to select a range of elements in matrix p0(i,j), i started from (at-9) to (at+9) and j started from 1 to n, (n=20) and the value of at is known. then do calculation for p0, but when i run it, it gives this error: "Attempted to access p0(NaN,2); index must be a positive integer or logical. Error in static (line 27) p0(i,j)=pin;" Please, I need to fix this error. thanks in advace.
clear all;clc;m=360; n=20; d=240; l=240;error=1e-4;RF=0.95;erin=0;error_alt=1*pi/180;deltheta=2*pi/m; delz=1/n;pin=1;s=(d*d)/(l*l);R=1/deltheta^2;cor=(deltheta^2/delz^2)*s;p0=zeros(m+1,n+1);p=zeros(m+1,n+1);p1=zeros(m+1,n+1);h=zeros(m+1);theta=zeros(m+1);ees=0.5; iter=1000; for i=1:m+1 theta(i)=(i-1)*deltheta; h(i)=1+ees*cos(theta(i));endsum(1)=0;alt1=atan(sqrt(1.0-ees^2)/ees);alt1=m+1-(alt1*180/pi);at=round(alt1);disp(at)for k=1:iter sumij=0; for j=2:n for i=(at-9):(at+9); p0(i,j)=pin; end end for j=2:n p0(m+1,j)=p0(1,j); end for j=2:n for i=1:m if (i==1) h1=1+ees*cos(theta(i)+0.5*deltheta); h2=1+ees*cos(theta(i)-0.5*deltheta); conh=h1-h2; const=2*(1+cor)*R*h(i)^3; A=3*h(i)^2*R*conh/const; C=h(i)^3*R/const; E=h(i)^3*R*cor/const; G=-3*sqrt(R)*conh/const; p(i,j)=(C+A/2)*p0(i+1,j)+(C-A/2)*p0(m,j)+E*(p0(i,j+1)+p0(i,j-1))+G; else conh=(h(i+1)-h(i-1))/2; const=2*(1+cor)*R*h(i)^3; A=3*h(i)^2*R*conh/const; C=h(i)^3*R/const; E=h(i)^3*R*cor/const; G=-3*sqrt(R)*conh/const; p(i,j)=(C+A/2)*p0(i+1,j)+(C-A/2)*p0(i-1,j)+E*(p0(i,j+1)+p0(i,j-1))+G; end if p(i,j)<0 p(i,j)=0.0; sumij=sumij+p(i,j); else end end end for j=2:n p(m+1,j)=p(1,j); end for j=2:n for i=(at-9):(at+9) p(i,j)=p0(i,j); end end sum(k+1)=sumij; percentage=abs(sum(k+1)-sum(k))/abs(sum(k+1)); if percentage < 0.0001 for j=1:n for i=1:m p1(i,j)=p(i,j); if p1(i,j) < 0 p1(i,j)=0; end end end else for j=2:n for i=1:m p0(i,j)=RF*p(i,j)+(1-RF)*p0(i,j); end end end wr=0; for j=1:n+1 wi=0.0; for i=1:m wi=wi+p1(i,j)*cos(theta(i))*deltheta; end wr=wr+wi*delz; end wt=0.0; for j=1:n+1 wi=0.0; for i=1:m wi=wi+p1(i,j)*sin(theta(i))*deltheta; end wt=wt+wi*delz; end w=sqrt(wr^2+wt^2); alt=atan(wt/wr); alt=m+1+(alt*180/pi); erin=abs(alt-alt1)/abs(alt); if erin<error_alt kjk=fprintf(1,'%4i\r\n',alt); break else alt1=0.5*alt+(1-0.5)*alt1; at=round(alt1); endendy=k;figure(1)surf(p1)figure(2)plot(p)
Best Answer