function new_crappy_vector = timestep2(num_points, x, y, Phi, h);

%Define wave function in polar coords
r=exp(y);
theta=-x;

%Initialize lots of derivatives
dr_dj=zeros(1, num_points);
dtheta_dj=zeros(1, num_points);
ds_dj=zeros(1, num_points);
dPhi_dj=zeros(1, num_points);
dr_ds=zeros(1, num_points);
dtheta_ds=zeros(1, num_points);
dPhi_ds=zeros(1, num_points);

% To make the derivatives work continuously along the circle,
%I used x = mod(x-1, num_points) + 1, so that if x is greater
%than num_points, it will curl back on itself. That is good. 

%Define the derivatives
for i=1:num_points
	dr_dj(i)=r(mod(i+1-1, num_points)+1)-r(i)  +  r(i)-r(mod(i-1-1, num_points)+1);
	dtheta_dj(i)=theta(mod(i+1-1, num_points)+1)-theta(i)  +  theta(i)-theta(mod(i-1-1, num_points)+1);
	dPhi_dj(i)=Phi(mod(i+1-1, num_points)+1)-Phi(i)  +  Phi(i)-Phi(mod(i-1-1, num_points)+1);
end

for i=1:num_points
    ds_dj(i)=sqrt((dr_dj(i)).^2+(dtheta_dj(i)).^2);
    dr_ds(i)=dr_dj(i)./ds_dj(i);
    dtheta_ds(i)=dtheta_dj(i)./ds_dj(i);
    dPhi_ds(i)=dPhi_dj(i)./ds_dj(i);
end

%calculate the dPhi_dn nasty integral

R=zeros(num_points, num_points);
alpha=zeros(num_points, num_points);
s=zeros(num_points, num_points);

for p1=1:num_points
for p2=1:num_points
	R(p1, p2)=sqrt((r(p1))^2+(r(p2))^2-2*r(p1)*r(p2)*cos(theta(p1)-theta(p2)));
	
	if p2 > p1
		for m = p1:p2-1
            s(p1, p2)=s(p1, p2)+(theta(m)-(theta(mod(m+1-1, num_points)+1)))*(r(m)+r(m+1))/2;
		end
	elseif p1 > p2
		for m = 1:p2-1
            s(p1, p2)=s(p1, p2)+(theta(m)-(theta(mod(m+1-1, num_points)+1)))*(r(m)+r(m+1))/2;
		end
		for m = p1: num_points-1
            s(p1, p2)=s(p1, p2)+(theta(m)-(theta(mod(m+1-1, num_points)+1)))*(r(m)+r(m+1))/2;
		end
		s(p1, p2)=s(p1, p2)+(theta(num_points)-(theta(1)-2*pi))*(r(num_points)+r(1))/2;
    else
        s(p1, p2) = 0;
    end

    s_N=s(1, num_points) + s(num_points, 1);

    if (p1 > p2) | (p1 < p2)
        alpha(p1, p2)=atan((r(p2)*sin(theta(p2))-r(p1)*sin(theta(p1)))/(r(p2)*cos(theta(p2))-r(p1)*cos(theta(p1))));
    else
        alpha(p1, p2)=atan((dtheta_dj(p1)*r(p1)*sin(theta(p1))+dr_dj(p1)*sin(theta(p1)))/(dr_dj(p1)*cos(theta(p1))-dtheta_dj(p1)*r(p1)*sin(theta(p1))));
    end
end
end


A=zeros(num_points, num_points);

for p1=1:num_points
	
    %Define A for the equation A*dPhi_dn=B
	for k=p1+1:p1+6
		p2 = mod(k-1, num_points)+1;
		A(p1,p2) = log(R(p1, p2)/s(p1, p2))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
		A(p1,p2) = A(p1, p2)+log(s(p1, p2))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
	end
	for k=p1+7:num_points+p1-7
        p2 = mod(k-1, num_points)+1;
        A(p1, p2) = log(R(p1, p2))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
	end
	for k=num_points+p1-6:num_points+p1-2
        p2 = mod(k-1, num_points)+1;
        A(p1, p2) = log(R(p1, p2)/(s_N-s(p1, p2)))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
        A(p1, p2) = A(p1, p2)+log(s_N-s(p1, p2))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
	end
    p2 = mod(p1-1-1, num_points)+1;
    A(p1, p2) = log(R(p1, p2)/(s_N-s(p1, p2)))*(s(p1, mod(p2+1-1, num_points)+1)-s(p1, p2));
    A(p1, p2) = A(p1, p2)+log(s_N-s(p1, p2))*(s_N-s(p1, p2));

end

% Solve the matrix equation to get dPhi_dn(i)

B = zeros(num_points,1);
for i=1:num_points
    for j=1:num_points
        B(i)=Phi(j)*(alpha(i, mod(j+1-1, num_points)+1)-alpha(i, j));
    end
end
A_inverse=inv(A);

dPhi_dn=zeros(1, num_points);
dPhi_dn = A_inverse*(B);
dPhi_dn=dPhi_dn'-2*pi*Phi;

%Finally do some DE calculation
%2nd Order Runge-Kutta. 

p_s = 0.146;

new_r_1 = r + (h/2)*((r.^2).*dr_ds.*dPhi_ds + (r.^2).*dtheta_ds.*dPhi_dn);
new_theta_1 = theta + (h/2)*(r.^2.*dtheta_ds.*dPhi_ds - r.*dr_ds.*dPhi_dn);
new_Phi_1 = Phi + (h/2)*(-p_s - log(r) + (1/2)*(r.^2).*((dPhi_ds).^2 + (dPhi_dn).^2));

new_r_2 = r + (h)*((new_r_1.^2).*dr_ds.*dPhi_ds + (new_r_1.^2).*dtheta_ds.*dPhi_dn);
new_theta_2 = theta + (h)*(new_r_1.^2.*dtheta_ds.*dPhi_ds - new_r_1.*dr_ds.*dPhi_dn);
new_Phi_2 = Phi + (h)*(-p_s - log(new_r_1) + (1/2)*(new_r_1.^2).*((dPhi_ds).^2 + (dPhi_dn).^2));

new_x = -new_theta_2;
new_y = log(new_r_2);
new_Phi = new_Phi_2;

new_crappy_vector = [new_x, new_y, new_Phi];

save('h:\scientific computing\project\xy_wave_2', 'num_points', 'new_x', 'new_y', 'new_Phi')

%subplot(2, 1, 1); plot(x, y);
%subplot(2, 1, 2); plot(new_x, new_y);