%Paramaterized boundaries
%By Brad Nelson 2/8/12
%m: number of equispaced periodic trapezoid nodes to use

function [z,dzdx,n,k]=param(m)

% parameters for circle
% r=@(x) 1; % 2pi periodic function, to be supplied by user
% drdx=@(x) 0; % above function's derivative.
% d2rdx2=@(x) 0; 

%parameters for problem
r=@(x) 1+0.3*cos(3*x); % 2pi periodic function, to be supplied by user
drdx=@(x) -0.9*sin(3*x); % above function's derivative.
d2rdx2=@(x) -2.7*cos(3*x); % 2nd deriv

%polar to cartesian coordinates
zfn=@(s) [cos(s).*r(s) ; sin(s).*r(s)];
dzdxfn=@(s) [-sin(s).*r(s)+cos(s).*drdx(s); cos(s).*r(s)+sin(s).*drdx(s)];
d2zd2xfn=@(s) [-cos(s).*r(s)-2*sin(s).*drdx(s)+cos(s).*d2rdx2(s);...
    -sin(s).*r(s)+2*cos(s).*drdx(s)+sin(s).*d2rdx2(s)];

s=0:2*pi/m:2*pi*(1-1/m); %sets up equispaced nodes on [0,2pi)

z=zfn(s); %position
dzdxvec=dzdxfn(s); %derivative at position
n=zeros(size(dzdxvec));

%creates vector orthogonal to z
n(1,:)=dzdxvec(2,:); 
n(2,:)=-dzdxvec(1,:);

for i=1:numel(n(1,:))
n(:,i)=n(:,i)./norm(n(:,i)); %normalizes n
end

dzdx=zeros(1,numel(dzdxvec(1,:)));

for i=1:numel(dzdxvec(1,:));
dzdx(i)=norm(dzdxvec(:,i)); %takes magnitude of derivative
end

k=-sum(n.*d2zd2xfn(s))./dzdx.^2; %curvature

% scatter3(z(1,:),z(2,:),k) %can use to test curvature for shape
% hold on
% scatter3(z(1,:),z(2,:),zeros(1,numel(z(1,:))),'x')  %actual 2d shape
% print -depsc 4a.eps;