/* file cylinder.mac 
Edwin L Woollett, 2008,
    woollett@charter.net
    http://www.csulb.edu/~woollett
    See tutorial notes:
     Maxima by Example, Ch. 6: Differential Calculus for
     a discussion of this method.
    
    this file sets up and uses
    cylindrical coordinates: (rho, phi, z) --> (rh, p, z).
    
  batch(sphere)  sets up environment and calculates
    the laplacian(f),  divergence(Bvec), and the
     (rho, phi, z) components of grad(f) and curl(Bvec) 
*/

"  cylindrical coordinates (rho, phi, z ) = (rh, p, z)  "$
   
" replacement rules x,y,z to rh, p, z  "$

c3rule : [x = rh*cos(p), y = rh*sin(p)  ]$
c3sub(expr) := (subst(c3rule,expr),trigsimp(%%) )$

rhxy : sqrt(x^2 + y^2)$

assume(rh > 0)$

" partial derivatives of rho and phi wrt x and y "$

drhdx : (diff(rhxy,x),c3sub(%%) );
drhdy : (diff(rhxy,y),c3sub(%%) );

dpdx : (diff( atan(y/x) ,x),c3sub(%%) );
dpdy : (diff( atan(y/x) ,y),c3sub(%%) );

" tell Maxima rh=rh(x,y) and p = p(x,y) "$


( gradef(rh,x,drhdx),gradef(rh,y,drhdy) )$
( gradef(p,x,dpdx),gradef(p,y,dpdy) )$

" tell Maxima to treat scalar function f as an "$
"    explicit function of (rh,p,z) "$

depends(f,[rh,p,z]);

"  calculate the Laplacian of the scalar function f(rh,p,z) "$
"  using the cartesian definition "$

( diff(f,x,2) + diff(f,y,2) + diff(f,z,2),trigsimp(%%),
          multthru(%%) );
                 
grind(%)$                  

" prepare to find derivatives which involve vectors. "$


"  cross product rule when using lists for vectors "$

lcross(u,v) := ( 
         ( u[2]*v[3] - u[3]*v[2] )*xu +
          ( u[3]*v[1] - u[1]*v[3]  )*yu +
          ( u[1]*v[2] - u[2]*v[1] )*zu
          
          )$
          
" cross product of parallel vectors is zero "$
              
lcross([a,b,c],[n*a,n*b,n*c]);
          
          
" a function we can use with map "$
          
apcr(ll) := apply('lcross, ll)$


" 3d cartesian unit vectors using lists "$

( xu : [1,0,0], yu : [0,1,0], zu : [0,0,1] )$   


" orthonormality checks on cartesian unit vectors "$

[ xu.xu, yu.yu, zu.zu, xu.yu, xu.zu, yu.zu ];


" low tech check of cross products of cartesian unit vectors"$

  lcross(xu,yu) - zu;

  lcross(yu,zu) - xu;

  lcross(zu,xu) - yu;
  
  [lcross(xu,xu),lcross(yu,yu),lcross(zu,zu)];
  
  
" high tech checks of cross products of cartesian unit vectors"$

  map('apcr,[[xu,yu],[yu,zu],[zu,xu]]) - [zu,xu,yu];
  
  map('apcr,[[xu,xu],[yu,yu],[zu,zu]]) ;  
  

" cylindrical coordinate unit vectors rho-hat, phi-hat  "$

rhu : cos(p)*xu + sin(p)*yu;
pu : -sin(p)*xu + cos(p)*yu;

" orthonormality checks on unit vectors "$

([ rhu.rhu, pu.pu, zu.zu, rhu.pu, rhu.zu, pu.zu], trigsimp(%%) );

" low tech check of cross products "$

 ( lcross(rhu,pu), trigsimp(%%) ) - zu;

 ( lcross(pu,zu), trigsimp(%%) ) - rhu;

 ( lcross(zu,rhu), trigsimp(%%) ) - pu;

[lcross(rhu,rhu), lcross(pu,pu)];

" high tech checks of cross products "$

  ( map('apcr,[[ rhu,pu], [pu,zu], [zu,rhu] ] ), trigsimp(%%) ) - [zu,rhu,pu];
  
   map('apcr, [ [rhu,rhu], [pu,pu] ] ) ;


" cartesian def. of gradient of scalar f "$

fgradient : diff(f,x)*xu + diff(f,y)*yu + diff(f,z)*zu $

" rho, phi, and z components of grad(f) "$

fgradient_rh : (rhu . fgradient, trigsimp(%%) );
fgradient_p : ( pu . fgradient, trigsimp(%%) );
fgradient_z : ( zu . fgradient, trigsimp(%%) );

" divergence of bvec "$

bvec : bx*xu + by*yu + bz*zu;

" two equations which relate cylindrical components"$
"   of bvec to the cartesian components "$

eq1 : brh = rhu.bvec;
eq2 : bp = pu.bvec;

" invert these equations "$

sol : (linsolve([eq1,eq2],[bx,by ]), trigsimp(%%) );

[bx,by] : map('rhs,sol)$

" tell Maxima to treat cylindrical components as "$
"   explicit functions of (rh,p,z) "$

depends([brh,bp,bz],[rh,p,z]);

" calculate the divergence of bvec "$

bdivergence : ( diff(bx, x) + diff(by, y) + diff(bz, z),trigsimp(%%),
                         multthru(%%) );

" cartesian definition of curl(vec)   "$

bcurl : (diff(bz,y) - diff(by,z) )*xu + (diff(bx,z) - diff(bz,x) )*yu +
                 (diff(by,x) - diff(bx,y) )*zu$
                 
" find cylindrical components of curl(bvec) "$

bcurl_rh : (rhu.bcurl,trigsimp(%%),multthru(%%) );
bcurl_p : (pu.bcurl,trigsimp(%%),multthru(%%) );
bcurl_z : (zu.bcurl,trigsimp(%%),multthru(%%) );



/*
(%i1) batch(cylinder);
batching #pc:/work2/cylinder.mac
(%i2)           cylindrical coordinates (rho,phi,z ) = (rh,p,z)  
(%i3)                 replacement rules x,y,z to rh,p,z  
(%i4)               c3rule : [x = rh cos(p), y = rh sin(p)]
(%i5)         c3sub(expr) := (subst(c3rule, expr), trigsimp(%%))
                                          2    2
(%i6)                        rhxy : sqrt(y  + x )
(%i7)                           assume(rh > 0)
(%i8)           partial derivatives of rho and phi wrt x and y 
(%i9)                 drhdx : (diff(rhxy, x), c3sub(%%))
(%o9)                               cos(p)
(%i10)                drhdy : (diff(rhxy, y), c3sub(%%))
(%o10)                              sin(p)
                                       y
(%i11)               dpdx : (diff(atan(-), x), c3sub(%%))
                                       x
                                     sin(p)
(%o11)                             - ------
                                       rh
                                       y
(%i12)               dpdy : (diff(atan(-), y), c3sub(%%))
                                       x
                                    cos(p)
(%o12)                              ------
                                      rh
(%i13)               tell Maxima rh=rh(x,y) and p = p(x,y) 
(%i14)           (gradef(rh, x, drhdx), gradef(rh, y, drhdy))
(%i15)             (gradef(p, x, dpdx), gradef(p, y, dpdy))
(%i16)           tell Maxima to treat scalar function f as an 
(%i17)                    explicit function of (rh,p,z) 
(%i18)                      depends(f, [rh, p, z])
(%o18)                           [f(rh, p, z)]
(%i19)      calculate the Laplacian of the scalar function f(rh,p,z) 
(%i20)                   using the cartesian definition 
(%i21) (diff(f, z, 2) + diff(f, y, 2) + diff(f, x, 2), trigsimp(%%), 
                                                                  multthru(%%))
                                   2
                                  d f
                            df    ---
                            ---     2    2     2
                            drh   dp    d f   d f
(%o21)                      --- + --- + --- + ----
                            rh      2     2      2
                                  rh    dz    drh
(%i22)                             grind(%)
'diff(f,rh,1)/rh+'diff(f,p,2)/rh^2+'diff(f,z,2)+'diff(f,rh,2)$
(%i23)        prepare to find derivatives which involve vectors. 
(%i24)          cross product rule when using lists for vectors 
(%i25) lcross(u, v) := (u  v  - u  v ) zu + (u  v  - u  v ) yu
                         1  2    2  1         3  1    1  3
                                                           + (u  v  - u  v ) xu
                                                               2  3    3  2
(%i26)             cross product of parallel vectors is zero 
(%i27)                lcross([a, b, c], [n a, n b, n c])
(%o27)                                 0
(%i28)                  a function we can use with map 
(%i29)                  apcr(ll) := apply('lcross, ll)
(%i30)               3d cartesian unit vectors using lists 
(%i31)         (xu : [1, 0, 0], yu : [0, 1, 0], zu : [0, 0, 1])
(%i32)          orthonormality checks on cartesian unit vectors 
(%i33)      [xu . xu, yu . yu, zu . zu, xu . yu, xu . zu, yu . zu]
(%o33)                        [1, 1, 1, 0, 0, 0]
(%i34)     low tech check of cross products of cartesian unit vectors
(%i35)                        lcross(xu, yu) - zu
(%o35)                             [0, 0, 0]
(%i36)                        lcross(yu, zu) - xu
(%o36)                             [0, 0, 0]
(%i37)                        lcross(zu, xu) - yu
(%o37)                             [0, 0, 0]
(%i38)         [lcross(xu, xu), lcross(yu, yu), lcross(zu, zu)]
(%o38)                 [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
(%i39)    high tech checks of cross products of cartesian unit vectors
(%i40)     map('apcr, [[xu, yu], [yu, zu], [zu, xu]]) - [zu, xu, yu]
(%o40)                 [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
(%i41)            map('apcr, [[xu, xu], [yu, yu], [zu, zu]])
(%o41)                 [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
(%i42)       cylindrical coordinate unit vectors rho-hat, phi-hat  
(%i43)                    rhu : sin(p) yu + cos(p) xu
(%o43)                        [cos(p), sin(p), 0]
(%i44)                    pu : cos(p) yu - sin(p) xu
(%o44)                       [- sin(p), cos(p), 0]
(%i45)               orthonormality checks on unit vectors 
(%i46) ([rhu . rhu, pu . pu, zu . zu, rhu . pu, rhu . zu, pu . zu], 
                                                                  trigsimp(%%))
(%o46)                        [1, 1, 1, 0, 0, 0]
(%i47)                 low tech check of cross products 
(%i48)               (lcross(rhu, pu), trigsimp(%%)) - zu
(%o48)                             [0, 0, 0]
(%i49)               (lcross(pu, zu), trigsimp(%%)) - rhu
(%o49)                             [0, 0, 0]
(%i50)               (lcross(zu, rhu), trigsimp(%%)) - pu
(%o50)                             [0, 0, 0]
(%i51)                [lcross(rhu, rhu), lcross(pu, pu)]
(%o51)                      [[0, 0, 0], [0, 0, 0]]
(%i52)                high tech checks of cross products 
(%i53) (map('apcr, [[rhu, pu], [pu, zu], [zu, rhu]]), trigsimp(%%))
                                                                - [zu, rhu, pu]
(%o53)                 [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
(%i54)                map('apcr, [[rhu, rhu], [pu, pu]])
(%o54)                      [[0, 0, 0], [0, 0, 0]]
(%i55)              cartesian def. of gradient of scalar f 
(%i56)     fgradient : diff(f, z) zu + diff(f, y) yu + diff(f, x) xu
(%i57)               rho, phi, and z components of grad(f) 
(%i58)          fgradient_rh : (rhu . fgradient, trigsimp(%%))
                                      df
(%o58)                                ---
                                      drh
(%i59)           fgradient_p : (pu . fgradient, trigsimp(%%))
                                      df
                                      --
                                      dp
(%o59)                                --
                                      rh
(%i60)           fgradient_z : (zu . fgradient, trigsimp(%%))
                                      df
(%o60)                                --
                                      dz
(%i61)                        divergence of bvec 
(%i62)                   bvec : bz zu + by yu + bx xu
(%o62)                           [bx, by, bz]
(%i63)         two equations which relate cylindrical components
(%i64)                 of bvec to the cartesian components 
(%i65)                      eq1 : brh = rhu . bvec
(%o65)                    brh = by sin(p) + bx cos(p)
(%i66)                       eq2 : bp = pu . bvec
(%o66)                    bp = by cos(p) - bx sin(p)
(%i67)                      invert these equations 
(%i68)       sol : (linsolve([eq1, eq2], [bx, by]), trigsimp(%%))
(%o68)    [bx = brh cos(p) - bp sin(p), by = brh sin(p) + bp cos(p)]
(%i69)                     [bx, by] : map('rhs, sol)
(%i70)          tell Maxima to treat cylindrical components as 
(%i71)                   explicit functions of (rh,p,z) 
(%i72)                depends([brh, bp, bz], [rh, p, z])
(%o72)            [brh(rh, p, z), bp(rh, p, z), bz(rh, p, z)]
(%i73)                 calculate the divergence of bvec 
(%i74) bdivergence : (diff(bz, z) + diff(by, y) + diff(bx, x), trigsimp(%%), 
                                                                  multthru(%%))
                                  dbp
                                  ---
                            brh   dp    dbz   dbrh
(%o74)                      --- + --- + --- + ----
                            rh    rh    dz    drh
(%i75)                cartesian definition of curl(vec)   
(%i76) bcurl : (diff(by, x) - diff(bx, y)) zu + (diff(bx, z) - diff(bz, x)) yu
                                               + (diff(bz, y) - diff(by, z)) xu
(%i77)             find cylindrical components of curl(bvec) 
(%i78)       bcurl_rh : (rhu . bcurl, trigsimp(%%), multthru(%%))
                                   dbz
                                   ---
                                   dp    dbp
(%o78)                             --- - ---
                                   rh    dz
(%i79)        bcurl_p : (pu . bcurl, trigsimp(%%), multthru(%%))
                                  dbrh   dbz
(%o79)                            ---- - ---
                                   dz    drh
(%i80)        bcurl_z : (zu . bcurl, trigsimp(%%), multthru(%%))
                                 dbrh
                                 ----
                                  dp    bp   dbp
(%o80)                         - ---- + -- + ---
                                  rh    rh   drh
(%o81)                       c:/work2/cylinder.mac
(%i82) 

*/