%%MATLAB ASSIGNMENT 2: DERIVATIVES REVISITED

%Welcome to MATLAB assignment 2, where calculus is made real

%The "real" numbers as you know them are an illusion, or a very useful 
%convention. The physical universe doesn't allow for arbitrary small
%values. So, for this assignment, we just agree that the smallest value
%we'll deal with is...

dx = 0.0001;

%..and use it to define everything else. For example, the differential:

d = @(f,x) f(x+dx) - f(x);

%d(f,x) computes df(x). You have seen the notation "du" in your 
%substituiton rule, but it was very loosely defined. Here, df(x) is a very
%concrete difference. 

%What is this difference good for? Well, the derivative is just the ratio
%df/dx. See for yourself:

f = @(x) x.^2; %we define a function f = x^2
should_be_six = d(f,3)/dx %this computes df/dx(3) -- what should the value be?

%This lab will revisit the derivative and its connection to the integral.

%%
%PART 1: THE FUNDAMENTAL THEOREM OF CALCULUS
%-------------------------------------------------------------------------
%Define a function int(f,a,b) that computes the integral numerically with
%step size equal to dx
%int(f,a,b) should return the integral of f from a to b
%HINT: use your knowledge from previous lab. You don't need to bother with
%left / right sums since the value of dx is small; an extra summand won't
%affect the precision much.

%You can put an anonymous function definition here, or make a function file.

%If you choose to do it with an anonymous function, put it here:
%int = @(f,a,b)  <---- uncomment and write the definition


%%                        
%Define a function f  = x^3. Compute its integral from 0 to 1 using your
%int() function. Verify that the value is what it should be.

%%
%Integral of f from 0 to 1 using your int():

%%
%True value (computed by hand):

%%
%Here, I define a new function in terms of yours:
df = @(x) d(f,x);

%%
%Compute the integral of df(x)/dx from 0 to 1, 2, 3 using your integral
%function. Verify that the values are what they should be according to FTC.
%   HINT: define a new function g = df(x)/dx, and integrate it.
%Alternatively, compute the integral of df, and divide by dx
%(MATLAB is stupid, and df./dx is not a valid expression in MATLAB)


%%
%Integral of df/dx from 0 to 1:

%%
%Should be equal to this according to FTC:

%%
%Integral of df/dx from 0 to 2:

%%
%Should be equal to this according to FTC:

%%
%Integral of df/dx from 0 to 3:

%%
%Should be equal to this according to FTC:


%%
%Let's go the other way around! 
%Define F = integral from 0 to x of f using your int function 
%(F is a function of x)
F = @(x) int(f,0,x);

%%
%Now compute d/dx(F(x)) for x = 1, 2, 3 USING the d(f,x) function defined 
%above, Verify that the values are what they should be.
%   NOTE: use the d() function to compute this numerically

%%
%d/dx(F(x)) at x = 1 is...


%%
%Should be equal to this according to FTC:


%%
%d/dx(F(x)) at x = 2 is...


%%
%Should be equal to this according to FTC:


%%
%d/dx(F(x)) at x = 3 is...


%%
%Should be equal to this according to FTC:


%%
%PART 2: THE SUBSTITUTION RULE
%-------------------------------------------------------------------------

%define a function u = x^2



%%
%Compute the integral of cos(u(x)) du/dx from 0 to sqrt(pi/6) using your the 
%d() function defined above and your int() function
%   HINT: define a function g  = cos(u(x)) du(x)/dx, (use the d() function!),
%   and plug it into your int().
%   This will compute integral of cos(u(x)) du/dx dx = cos(u) du
%   NOTE: don't differentiate / integrate by hand. 



%%
%Verify that the integral is what it should be according to the substitution rule.
%What should the value of the integral be? (now, compute it by hand):
%%
%
% $$\int_0^{\sqrt{\frac{\pi}{6}}} \cos(u(x)) \frac{du}{dx} dx = $$
%


%%
%PART 3: WE DIDN'T LIE, DERIVATIVES REALLY GIVE YOU TANGENTS
%-------------------------------------------------------------------------
figure('Name', 'Tangent and normal');

%Define  f = sqrt(x) and plot it between 0 and 4

hold on;
axis equal;      %sets axis scaling to 1:1 (so the figure is not stretched)

%%
%Use the segment command to draw the TANGENT line to the plot at x = 1.5
%   HINT: the slope is df/dx, so the direction vector is (dx, df)

%%
%Use the segment command to draw the NORMAL line to the plot at x = 1.5
%   HINT: the slope is -dx/df

%%
%PART 4: ...AND TANGENTS GIVE YOU DERIVATIVES
%-------------------------------------------------------------------------
%
%Here we illustrate that the slope of the secant lines approaches the
%derivative.
%
%%
%Set f(x) = sin(x) and plot it between 0 and pi

figure('Name', 'Derivative as limit of tangents');

%(put plot code here)

hold on;
axis equal;

%%
%Now, we'll plot N secant lines with a FOR loop
N = 20;

%The secant lines go through (x0, f(x0)) and (x1, f(x1))
%x0 will be fixed:
x0 = pi/3;
%x1 will take different values between MAX=2 and x0 in the loop
MAX = 2;
for i = 0:N-1 %i runs from 0 to N-1
    x1 =  MAX - (MAX-x0)*i/N;  
    %Compute the direction vector:
    %DX is the difference in X values:
    DX = x1 - x0;
    %DY is the difference in Y values:
    %DY =  %<----uncomment and put the value here
    %Use the segment command to plot the tangent line:
    
end

%%
%Output the last secant's slope DY/DX (this is last value it has taken in the loop) 
%and compare it to the tangent slope, that is, df/dx at x0 [uncomment lines below]

%last_secant_slope =  
%tangent_slope = 

%%
%PART 5: ENVELOPE
%-------------------------------------------------------------------------
figure('Name', 'Envelope');

%define f = x^2 and plot it between -1 and 1

%(put the plot code here)

hold on ;
axis equal;

%%
%Draw tangents at 100 points between -1 and 1
%   HINT: set up an array, and use a for loop to loop through it
%   NOTE: pick a nice value for L in the segment command - somethine
%   between 2 and 5, for example - whatever you think looks good



%We say that the parabola is an envelope of its tangent lines.


%%
%PART 6: PARAMETRIC CURVE AS AN ENVELOPE OF ITS TANGENTS
%-------------------------------------------------------------------------
dt = dx;
t = [0:dt:2*pi];

figure('Name', 'Parametric envelope');
%The following defines and plots an ellipse
x1 = @(t) 3*cos(t);
x2 = @(t) 2*sin(t);
plot(x1(t),x2(t));
axis equal;
hold on;

%%
%Draw tangents at 100 equally spaced values of t between 0 and 2*pi
%   HINT: set up an array, and use a for loop to loop through it;
%   a tangent line at t = t0 has direction vector [dx1(t0), dx2(t0) ]

%%
%PART 6: PARAMETRIC CURVE NORMALS
%-------------------------------------------------------------------------
figure('Name', 'Normals to a curve');
%The following defines and plots an ellipse
x1 = @(t) 2*cos(t);
x2 = @(t) 3*sin(t);
t=[0:0.01:2*pi];
plot(x1(t),x2(t));

axis equal;
hold on;

%%
%Draw normals at 100 equally spaced values of t between 0 and 2*pi
%   HINT: set up an array, and use a for loop to loop through it
%   NOTE: set the L argument of the segment command to be 5 or greater for
%   a nice picture

%%
%Notice the shape that is formed inside. Is it smooth? Describe it here:



%%PART 7: EVOLUTES ------> EXTRA CREDIT, WORTH 1/2 LAB
%This section is a study of the curve in the previous section.
%If you think about it, what your eyes see is the curve formed by
%intersection points of consecutive normals. That is, the evolute is the
%<b>envelope</b> of the normals. Now, we just have to do the opposite of
%what we did in parts 4 and 5: given a family of lines, we construct
%a curve to which they are all tangent.

%%
%The tedious part: given two lines:
%L0 goes through (x0, y0) and has direction vector (a0, b0)
%L1 goes through (x1, y1) and has direction vector (a1, b1)
%What are the coordinates of the intersection point P?
%P satisfies: P = [x0 y0] + t0[a0 b0] = [x1 y1] + t1[a1 b1] 
%for some constants t0, t1. Set up and solve the system of linear equations
%for t0 and t1; obtain P as a function of all other variables



%%
%Use the function P to define a function PN which, given functions 
%x(t) and y(t) and values t0, t1, computes the intersection of normals
%at t = t0 and t = t1
%Make the output of the function an array of size 2 containing the x and y
%coordinates of the intersection



%%
%Now we can use the function PN to compute the points we want to plot.
%(That's  the fun part! :) )
%First, we define a curve:
figure('Name', 'Evolute of an ellipse');
x = @(t) 2*cos(t);
y = @(t) 3*sin(t);
N=100;

%Define an array t of size N with numbers going from 0 to 2*pi
t = linspace(0,2*pi,N);
  
%Plot (x(t), y(t))
plot(x(t),y(t));
axis equal;
hold on;

ex=zeros(N-1,1); %allocate an array for evolute x-points
ey=zeros(N-1,1); %allocate an array for evolute x-points
 for i=1:N-1
        %compute the intersection of normal through t(i) and t(i+1)
        %intersection = (uncomment and put something here)
        %store the intersection point 
        %ex(i) = intersection(1);
        %ey(i) = intersection(2);
        %connect the intersection point to the point on the curve at t(i)
        %HINT: plot([a c], [b d]) connects (a, b) to (c, d)
        
 end
%plot the evolute
plot(ex, ey);

%%
%Extra Extra credit (worth 1/4 of the lab), plus you get to make something
%Wikipedia-article worthy

%%
%Some math background first.
%
%The method of construction of the evolute above reveals its interesting
%property: it is the locus of the centers of the osculating circles,
%the circles that "kiss" the curve (from Latin "osculare" - to kiss).
%
%OK, hold on, let's talk about what this means first.
%
%Just like there are tangent lines, there are tangent circles.
%Unlike a tangent line, a circle tangent to a point is not unque - you can 
%draw many, by changing the radius (line can be thought of as a tangent 
%circle with infinite radius). 
%
%If you draw a lot of circles and zoom at the common tangent point, you
%will see that the larger circles are going on the _outside_ of the curve,
%and the smaller circles are on the _inside_.
%The one circle that is in between in the osculating circle.
%It it the limit of the circles going through three points on the curve
%as they approach each other.
%
%Note for physics people: a car going around the curve with unit speed will 
%experience the same acceleration at a point (x, y) as the car going around
%the osculating circle at that point with the same speed.
%
%We say that the curve at that point and the circle have the same _curvature_.
%
%In other words, just as the line is the best linear approximation, the
%circle is the best quadratic approximation.
%---------------------------------------------------------------------------

%%
%The coding part starts here.
%
%Write a circ(cx, cy, r) function that plots a circle of radius r centered at (cx, cy) 
%HINT: see how the ellipse is plotted.
%

%
%Write a function dist(x0, y0, x1, y1) that computes the distance 
%from (x0, y0) to (x1, y1)
%

%   
%Time to get things moving!
F = figure('Name', 'Evolute animated!'); %store reference to the figure in variable F
for i = 1:N-1 %this sets up the animation loop. Each iteration will draw a new frame.
   %Plot the original curve (x(t), y(t)):
   %
   
   %
   axis equal; %plot formatting 
   axis ([-8 8 -6 6]);
   hold on;
   %Plot a part of the evolute from time t=t(1) to t = t(i) 
   %[NOTE: use the ex, ey arrays that are already computed
   %
   
   %
   %Compute the radius of curvature R at time t=i
   %(this is the distance from the point on the curve at time t=t(i) to the corresponding
   % point on the evolute); use your dist function
   %
   
   %
   %Draw a circle of radius R centered at the point on the evolute at time t(i):
   %[use your circ function]
   %
   
   %
   %Plot a segment connecting the center of the circle to the correspoing
   %point on the curve (to which it will be tangent):
   %
   
   %
   hold off; %allow the next frame to overwrite the figure
   pause(1/24); %pause before the next frame gets drawn

   %The code below saves the animation to an animated gif file.
   %Once you get the animation right, uncomment these lines
   %to get a nice .gif file that you can post on your Tumblr.
   %Submit the GIF file along with the assignment (ideally, edit the
   %published HTML file to include instead of the static figure).
   %------------------------------------------------------------------
   %[A,map] = rgb2ind(frame2im(getframe(F)),256);
   %if i == 1;
   %imwrite(A,map,'evolute.gif','gif','LoopCount',Inf,'DelayTime',1/24);
   %else
   %	imwrite(A,map,'evolute.gif','gif','WriteMode','append','DelayTime',1/24);
   %end      
end
   
