subroutine.m'ҼsmBINcfunction [X,Y,K,L,nmax,E1,E2,PE1,PE2,epsilon,l0] = subroutine(T, deltat, x0, y0, k0, c0, G, variablec,g) nmax = floor(T/deltat); l0=abs(G*x0*k0*(c0^2*4*pi^2*(x0^2+y0^2)^2-G^2*x0^2)^-.5); % for initially constant in y % =.45 for good spiral X(1)=x0; Y(1)=y0; K(1)=k0; L(1)=l0; for n=2:nmax kappa=sqrt((K(n-1))^2+(L(n-1))^2); r=sqrt((X(n-1))^2+(Y(n-1))^2); u=-G*Y(n-1)/(2*pi*r^2); v=G*X(n-1)/(2*pi*r^2); A=G/(2*pi*r^4); B=2*X(n-1)*Y(n-1); C=(Y(n-1))^2-(X(n-1))^2; %if variablec %U2 = (G/(2*pi*r))^2; %c = c0*sqrt(1-((g-1)*U2)/(2*c0^2)); %a = -((g-1)*U2*kappa)/(2*c*r^2); %else c = c0; a = 0; %end X(n)=deltat*(c*K(n-1)/kappa+u)+X(n-1); Y(n)=deltat*(c*L(n-1)/kappa+v)+Y(n-1); K(n)=deltat*(a*X(n-1) - A*(B*K(n-1)+C*L(n-1)))+K(n-1); L(n)=deltat*(a*Y(n-1) + A*(B*L(n-1)-C*K(n-1)))+L(n-1); E2(n)=c*kappa+u*K(n-1)+v*L(n-1); end E2(1)=E2(2); E1=K.*Y-L.*X; PE1=abs(100*(max(E1)-min(E1))/mean(E1)); PE2=abs(100*(max(E2)-min(E2))/mean(E2)); epsilon=abs(G)/(2*pi*c0*abs(y0));