kill(all)$
load(lrats)$
load(vect)$

Calculation of curvature (= 1/radius ) ,(x,y) and (r,θ), 'version 1.0a'
peter.vlasschaert@gmail.com,01/10/20018

rem :  s = r*θ (measured along the circle or the curve)

 curvature

with use Figure 2 
rem : K,K1 are points 

Definition : Curvature (see : Figure 2)
===============================

p1: K = abs('diff(φ,s));

We want to calculation : 'diff(φ,s,1)

p2:part(rhs(p1),1);

   ?  arclength, pythagoras 
=======================

-------------------------------------------------------------------------------------
| Δs|^2 = |Δy|^2+|Δx|^2   =>   limit  : ds = sqrt(dx^2+dy^2)
-------------------------------------------------------------------------------------
math notation ,   dy/dx = y ' = yprime
-------------------------------------------------------------------------------------
math notation ,   ds/dx = sqrt(1+y'^2)

p3:'diff(s,x) =sqrt(1+'diff(y,x)^2);

definition : tan(φ)  = diff(y,x)
=================================
-----------------------------------------------------------
notation :  y = f(x)
notation : y' = diff(y,x) = tan(φ) = yprime
-----------------------------------------------------------

depends([yprime],x);

p4: tan(φ) = yprime;

p4a:solve(p4,φ)[1];

differentation : yprime with respect to x
------------------------------------------------------------

p5: 'diff(lhs(p4a),x) = diff(rhs(p4a),x);

p5a:yprime= 'diff(y,x);

p5b:subst([p5a],p5);

calculation : dφ/ds;
================

p6:p2= 'diff(φ,x) /'diff(s,x);

use : 'p3' and 'p5b'

p6a :subst([p3,p5b],p6);

final result: curvature
==================

p7:subst([p6a],p1);

use: curvature at any point  of a curve where second derivative
        (=diff(y,x,2)) exists and is continuous,only positive values
        of the root should be taken.

calculation of curvature : 'polar coordinates'
=====================================

r = r(θ)

depends([r],θ);

depends([r],θ);

polar coordinates : (x,y) => (r,θ)
==========================

p8:x=r*cos(θ);
p8a:y=r*sin(θ);

differentation x,y to θ , first time with respect to θ
=========================================

p9:'diff(lhs(p8),θ) =diff(rhs(p8),θ);
p9a:'diff(lhs(p8a),θ) =diff(rhs(p8a),θ);

yprime = y' = diff(y,x)
-------------------------------

p10:'diff(y,x)  = lhs(p9a)/lhs(p9);

p10a :lhs(p10) =  rhs(p9a)/rhs(p9);

p10b:p10a^2;

p10c:expand(p10b);

p10d:trigsimp(p10c);

p10e:1+p10d;

p10f:ratsimp(p10e);

p10g:trigsimp(p10f);

p10h:p10g^(3/2);

differentation x,y to θ ,'p9,p9a' with respect to θ
=========================================

p11:'diff(lhs(p8),θ,2) =diff(rhs(p8),θ,2);
p11a:'diff(lhs(p8a),θ,2) =diff(rhs(p8a),θ,2);

rem : d/dx(dy/dx) = d/dθ(dy/dθ/dx/dθ)*dθ/dx

find :d/dx(dy/dx)
==============

p12:diff((rhs(p9a)/rhs(p9)),θ)$

p12a:ratsimp(p12);

p12b:trigsimp(p12a);

p13:'diff(θ,x)=1/rhs(p9);

p14:p12b*rhs(p13);

p15:expand(p14);

p15a:trigsimp(p15);

problem : p10h we get complex values , ' result problem'

other way find result  :  curve ' represented' parametrically '
=================================================

use :parametric, 't'
================

depends([α,β],t);

α ,β  functions of t depent.

p16: x= α;
p16a:y = β;

find : dy/dx = dy/dt/(dx/dt) = dβ/dt/dα/dt.
-----------------------------------------------------------

p17:'diff(y,x) = 'diff(lhs(p16a),t)/ 'diff(lhs(p16),t);

p17a: lhs(p17) = subst([p16,p16a],rhs(p17));

p17a1:part(rhs(p17a),1);
p17a2:part(rhs(p17a),2);

find : d/dx(dy/dx) = d/dt(dy/dt/dx/dt)*dt/dx
------------------------------------------------------------

first: d/dt((dy/dt/dx/dt))
===================

p18:diff(rhs(p17a),t);
p18a:ratsimp(p18);

second : dt/dx = 1/(dx/dt)
=====================

p18b:part(rhs(p17),1);
p18c:part(rhs(p17),2);

p18d: 'diff(t,x) = 1/p18c;

p18e:lhs(p18d) = subst([p16],rhs(p18d));

thrid : first*second
================

p19:'diff(y,x,2) = p18a*rhs(p18e);

to find curvature K  :
================

p20: p7;

subst. 'p17a,p19' into 'p20'

p21:subst([p17a,p19],p20);

p21a:abs(p17a2)=p17a2;

p22:subst([p21a],p21);

study :simplification

p23:('diff(β,t,1))^2/('diff(α,t,1))^2+1;

p24:ratsimp(p23);

p24a:p24^(3/2);

p24b:subst([p21a],p24a);

final :result parametric representation of the curvature
==============================================

p25a:p16;
p25b:p16a;

notation : x=α(t),y=β(t)
--------------------------

definition : K  ,with 't' = parametric value
==================================

p26:part(rhs(p22),1);
p26a:part(rhs(p22),2,1);
p26b:part(rhs(p22),2,2) =p24b;

final equation : curvature 'parametrically defined curve'
===============================================

p27:lhs(p22) = p26/(p26a*rhs(p26b));

try equation 'p27'  to find , curvature for 'polar coordinates'
=================================================

p28:subst([t=θ],p27);

polar coordinates :
================
  α => r*cos(θ)
  β => r*sin(θ)

p28a:subst([α = r*cos(θ),β = r*sin(θ)],p28);

remove : abs  values, why  ?
========================
  you get lot problems to simplify equation 'p28' or 'p28a'

p28b: abs(part(rhs(p28a),1)) = part(rhs(p28a),1,1);
p28c:subst([p28b],p28a);

excute : differentation :
----------------------------------

p28d:ev(p28c,diff);

final equation : curvature 'polar coordinates'
=====================================

p28e:trigsimp(p28d)$
p28f: lhs(p28e) = abs(part(rhs(p28e),1));

final equation : radius of curvature
==============================

definition :  R = 1/K

p29: R= 1/rhs(p28f);

-----------------------------------------------------------------------------------------------------------------------------------------------------------------