Steve Roper

180 Reputation

8 Badges

6 years, 136 days

MaplePrimes Activity


These are replies submitted by Steve Roper

Updated 26 January 2022 for compatibility with Maple 2021.2:

Antenna_arrays_2021.2.mw

Updated 26 January 2022 for compatibility with Maple 2021.2:

Friis_Transmission_Equation_2021.2.mw

@Samir Khan Thank you for your suggestion Samir,

Unfortunately any attempt to upload Antenna arrays.mw results in an error message.  The worksheet does take a while to process (on my MacBook), so perhaps it exceeds a timeout on the webserver?

@ecterrab Hi, Edgardo,

This version also includes the calculation of wave impedance and the Poynting Vector...

Examining Maxwell's EquationsNULL

Initialise

   

Maxwell's Equations

 

We define Maxwell's Equations in terms of vector functions in Cartesian coordinates (this is the simple case for the vacuum):

Maxwell–Faraday equation:

Maxwell_1 := Curl(E__field_(x, y, z, t)) = -(diff(B__flux_(x, y, z, t), t))

Physics:-Vectors:-Curl(E__field_(x, y, z, t)) = -(diff(B__flux_(x, y, z, t), t))

(2.1)

Ampère's circuital law (with Maxwell's addition):

Maxwell_2 := Curl(H__field_(x, y, z, t)) = diff(D__flux_(x, y, z, t), t)

Physics:-Vectors:-Curl(H__field_(x, y, z, t)) = diff(D__flux_(x, y, z, t), t)

(2.2)

Gauss' Law:

Maxwell_3 := Divergence(D__flux_(x, y, z, t)) = 0

Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)) = 0

(2.3)

Gauss' Law for Magnetism:

Maxwell_4 := Divergence(B__flux_(x, y, z, t)) = 0

Physics:-Vectors:-Divergence(B__flux_(x, y, z, t)) = 0

(2.4)

``

Constitutive Relations

 

Eq_1 := D__flux_(x, y, z, t) = epsilon*E__field_(x, y, z, t)

D__flux_(x, y, z, t) = varepsilon*E__field_(x, y, z, t)

(3.1)

Eq_2 := B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

(3.2)

``

General Solution for Maxwell's Equations

 

Finding a solution for the Electric Field

 

We need an expression for the Curl of H, so we start by taking the Curl of both sides of Maxwell_1:

Curl(Maxwell_1)

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(E__field_(x, y, z, t))) = -Physics:-Vectors:-Curl(diff(B__flux_(x, y, z, t), t))

(4.1.1)

Now substitute B for μH:

subs(Eq_2, Physics[Vectors][Curl](Physics[Vectors][Curl](E__field_(x, y, z, t))) = -Physics[Vectors][Curl](diff(B__flux_(x, y, z, t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(E__field_(x, y, z, t))) = -Physics:-Vectors:-Curl(diff(mu*H__field_(x, y, z, t), t))

(4.1.2)

Now substitute the Curl of H with an expression in D:

subs(diff(Maxwell_2, t), Physics[Vectors][Curl](Physics[Vectors][Curl](E__field_(x, y, z, t))) = -Physics[Vectors][Curl](diff(mu*H__field_(x, y, z, t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(E__field_(x, y, z, t))) = -mu*(diff(diff(D__flux_(x, y, z, t), t), t))

(4.1.3)

Now substitute E for D:

subs(isolate(Eq_1, E__field_(x, y, z, t)), Physics[Vectors][Curl](Physics[Vectors][Curl](E__field_(x, y, z, t))) = -mu*(diff(diff(D__flux_(x, y, z, t), t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(D__flux_(x, y, z, t)/varepsilon)) = -mu*(diff(diff(D__flux_(x, y, z, t), t), t))

(4.1.4)

Let's see what happens if we expand this:

expand(Physics[Vectors][Curl](Physics[Vectors][Curl](D__flux_(x, y, z, t)/varepsilon)) = -mu*(diff(diff(D__flux_(x, y, z, t), t), t)))

Physics:-diff(Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)), z)*_k/varepsilon+Physics:-diff(Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)), y)*_j/varepsilon+Physics:-diff(Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)), x)*_i/varepsilon-(diff(diff(D__flux_(x, y, z, t), y), y))/varepsilon-(diff(diff(D__flux_(x, y, z, t), x), x))/varepsilon-(diff(diff(D__flux_(x, y, z, t), z), z))/varepsilon = -mu*(diff(diff(D__flux_(x, y, z, t), t), t))

(4.1.5)

This is impressive!  Maple knows that V×(V×D) = V(V·D) - V2D

Now we can use Maxwell_3 to arrive at an expression in V2D because we know that V·D=0:

simplify(Physics[Vectors][Curl](Physics[Vectors][Curl](D__flux_(x, y, z, t)/varepsilon)) = -mu*(diff(diff(D__flux_(x, y, z, t), t), t)), {Maxwell_3})

(-(diff(diff(D__flux_(x, y, z, t), x), x))-(diff(diff(D__flux_(x, y, z, t), y), y))-(diff(diff(D__flux_(x, y, z, t), z), z)))/varepsilon = -mu*(diff(diff(D__flux_(x, y, z, t), t), t))

(4.1.6)

Now we have an expression that looks like a wave, we can substitute D for E:

subs(Eq_1, (-(diff(diff(D__flux_(x, y, z, t), x), x))-(diff(diff(D__flux_(x, y, z, t), y), y))-(diff(diff(D__flux_(x, y, z, t), z), z)))/varepsilon = -mu*(diff(diff(D__flux_(x, y, z, t), t), t)))

(-(diff(diff(varepsilon*E__field_(x, y, z, t), x), x))-(diff(diff(varepsilon*E__field_(x, y, z, t), y), y))-(diff(diff(varepsilon*E__field_(x, y, z, t), z), z)))/varepsilon = -mu*(diff(diff(varepsilon*E__field_(x, y, z, t), t), t))

(4.1.7)

simplify((-(diff(diff(varepsilon*E__field_(x, y, z, t), x), x))-(diff(diff(varepsilon*E__field_(x, y, z, t), y), y))-(diff(diff(varepsilon*E__field_(x, y, z, t), z), z)))/varepsilon = -mu*(diff(diff(varepsilon*E__field_(x, y, z, t), t), t)))

-(diff(diff(E__field_(x, y, z, t), x), x))-(diff(diff(E__field_(x, y, z, t), y), y))-(diff(diff(E__field_(x, y, z, t), z), z)) = -mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

(4.1.8)

Loose the "-" sign:

diff(diff(E__field_(x, y, z, t), x), x)+diff(diff(E__field_(x, y, z, t), y), y)+diff(diff(E__field_(x, y, z, t), z), z) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

diff(diff(E__field_(x, y, z, t), x), x)+diff(diff(E__field_(x, y, z, t), y), y)+diff(diff(E__field_(x, y, z, t), z), z) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

(4.1.9)

 

Finally, we can use the inert operator "%" to make the expression look more concise like this:

subs(Laplacian(E__field_(x, y, z, t)) = %Laplacian(E__field_(x, y, z, t)), diff(diff(E__field_(x, y, z, t), x), x)+diff(diff(E__field_(x, y, z, t), y), y)+diff(diff(E__field_(x, y, z, t), z), z) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t)))

%Laplacian(E__field_(x, y, z, t)) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

(4.1.10)

Finding a solution for the Magnetic Field

 

We need an expression for the Curl of E, so we start by taking the Curl of both sides of Maxwell_2:

Curl(Maxwell_2)

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(H__field_(x, y, z, t))) = Physics:-Vectors:-Curl(diff(D__flux_(x, y, z, t), t))

(4.2.1)

Now substitute D for ϵE:

subs(Eq_1, Physics[Vectors][Curl](Physics[Vectors][Curl](H__field_(x, y, z, t))) = Physics[Vectors][Curl](diff(D__flux_(x, y, z, t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(H__field_(x, y, z, t))) = Physics:-Vectors:-Curl(diff(varepsilon*E__field_(x, y, z, t), t))

(4.2.2)

Now substitute the Curl of E with an expression in B:

subs(diff(Maxwell_1, t), Physics[Vectors][Curl](Physics[Vectors][Curl](H__field_(x, y, z, t))) = Physics[Vectors][Curl](diff(varepsilon*E__field_(x, y, z, t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(H__field_(x, y, z, t))) = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t))

(4.2.3)

Now substitute H for B:

subs(isolate(Eq_2, H__field_(x, y, z, t)), Physics[Vectors][Curl](Physics[Vectors][Curl](H__field_(x, y, z, t))) = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t)))

Physics:-Vectors:-Curl(Physics:-Vectors:-Curl(B__flux_(x, y, z, t)/mu)) = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t))

(4.2.4)

expand(Physics[Vectors][Curl](Physics[Vectors][Curl](B__flux_(x, y, z, t)/mu)) = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t)))

Physics:-diff(Physics:-Vectors:-Divergence(B__flux_(x, y, z, t)), z)*_k/mu+Physics:-diff(Physics:-Vectors:-Divergence(B__flux_(x, y, z, t)), y)*_j/mu+Physics:-diff(Physics:-Vectors:-Divergence(B__flux_(x, y, z, t)), x)*_i/mu-(diff(diff(B__flux_(x, y, z, t), y), y))/mu-(diff(diff(B__flux_(x, y, z, t), x), x))/mu-(diff(diff(B__flux_(x, y, z, t), z), z))/mu = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t))

(4.2.5)

``

``

Now we can use Maxwell_4 to arrive at an expression in V2H because we know that V·B=0:

simplify(Physics[Vectors][Curl](Physics[Vectors][Curl](B__flux_(x, y, z, t)/mu)) = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t)), {Maxwell_4})

(-(diff(diff(B__flux_(x, y, z, t), x), x))-(diff(diff(B__flux_(x, y, z, t), y), y))-(diff(diff(B__flux_(x, y, z, t), z), z)))/mu = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t))

(4.2.6)

Now we have an expression that looks like a wave, we can substitute B for H:

subs(Eq_2, (-(diff(diff(B__flux_(x, y, z, t), x), x))-(diff(diff(B__flux_(x, y, z, t), y), y))-(diff(diff(B__flux_(x, y, z, t), z), z)))/mu = -varepsilon*(diff(diff(B__flux_(x, y, z, t), t), t)))

(-(diff(diff(mu*H__field_(x, y, z, t), x), x))-(diff(diff(mu*H__field_(x, y, z, t), y), y))-(diff(diff(mu*H__field_(x, y, z, t), z), z)))/mu = -varepsilon*(diff(diff(mu*H__field_(x, y, z, t), t), t))

(4.2.7)

simplify((-(diff(diff(mu*H__field_(x, y, z, t), x), x))-(diff(diff(mu*H__field_(x, y, z, t), y), y))-(diff(diff(mu*H__field_(x, y, z, t), z), z)))/mu = -varepsilon*(diff(diff(mu*H__field_(x, y, z, t), t), t)))

-(diff(diff(H__field_(x, y, z, t), x), x))-(diff(diff(H__field_(x, y, z, t), y), y))-(diff(diff(H__field_(x, y, z, t), z), z)) = -varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

(4.2.8)

Loose the "-" sign:

diff(diff(H__field_(x, y, z, t), x), x)+diff(diff(H__field_(x, y, z, t), y), y)+diff(diff(H__field_(x, y, z, t), z), z) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

diff(diff(H__field_(x, y, z, t), x), x)+diff(diff(H__field_(x, y, z, t), y), y)+diff(diff(H__field_(x, y, z, t), z), z) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

(4.2.9)

 

Finally, we can use the inert operator "%" to make the expression look more concise like this:

subs(Laplacian(H__field_(x, y, z, t)) = %Laplacian(H__field_(x, y, z, t)), diff(diff(H__field_(x, y, z, t), x), x)+diff(diff(H__field_(x, y, z, t), y), y)+diff(diff(H__field_(x, y, z, t), z), z) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t)))

%Laplacian(H__field_(x, y, z, t)) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

(4.2.10)

General PDE solutions

 

So the general solutions to Maxwell's Equations are:

E__fieldPDE := %Laplacian(E__field_(x, y, z, t)) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

%Laplacian(E__field_(x, y, z, t)) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

(4.3.1)

H__fieldPDE := %Laplacian(H__field_(x, y, z, t)) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

%Laplacian(H__field_(x, y, z, t)) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

(4.3.2)

Finding the equation for the magnitude of the Electric Field

 

infolevel[pdsolve] := 3

3

(4.4.1)

pdsolve(value(E__fieldPDE), E__field_(x, y, z, t))

PDESolStruc(E__field_(x, y, z, t) = _F1(x)*_F2(y)*_F3(z)*_F4(t), [{diff(diff(_F1(x), x), x) = _c[1]*_F1(x), diff(diff(_F2(y), y), y) = _c[2]*_F2(y), diff(diff(_F3(z), z), z) = _c[3]*_F3(z), diff(diff(_F4(t), t), t) = _F4(t)*_c[1]/(mu*varepsilon)+_F4(t)*_c[2]/(mu*varepsilon)+_c[3]*_F4(t)/(mu*varepsilon)}])

(4.4.2)

pdsolve has found four functions by the separation of variables x, y, z and t.  We can see that the double derivative of functions _F1(x), _F2(y) & _F3(z)  must equal the original function, multiplied by some constant.

The sin function satisfies this requirement:

diff(sin(x), x, x)

-sin(x)

(4.4.3)

And so does the exponential in 'i':

diff(exp(I*x), x, x)

-exp(I*x)

(4.4.4)

The benefit of using the exponential form will be that the product _F1(x),*_F2(y) *_F3(z)*_F4(t) becomes a single exponential with the exponent being a simple arithmetic sum, whereas a solution based on sin, would be a product with the form sin(x)*sin(y)*sin(z)*sin(t) and this would be difficult to simplify.

 

We will multiply the x, y and z coordinates by constants kx, ky and kz respectively and we will multiply the time coordinate by a constant ω:``

"`E__wave_`(x,y,z,t):=E*(e)^(i*(`k__x`*x))*(e)^(i*(`k__y`*y))*(e)^(i*(`k__z`*z))*(e)^(i*(-omega*t));"

proc (x, y, z, t) options operator, arrow, function_assign; Physics:-`*`(E, exp(Physics:-`*`(I, Physics:-`*`(k__x, x))), exp(Physics:-`*`(I, Physics:-`*`(k__y, y))), exp(Physics:-`*`(I, Physics:-`*`(k__z, z))), exp(Physics:-`*`(I, -Physics:-`*`(omega, t)))) end proc

(4.4.5)

Now test to see whether this wave equation is indeed a solution to the electric field PDE:

pdetest(E__field_(x, y, z, t) = E__wave*_, E__fieldPDE)

0

(4.4.6)

It worked!

 

Substitute the expression for the electric field wave into the electric field PDE:

subs(E__field_(x, y, z, t) = E__wave_(x, y, z, t), E__fieldPDE)

%Laplacian(E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t)) = mu*varepsilon*(diff(diff(E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t), t), t))

(4.4.7)

value(%Laplacian(E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t)) = mu*varepsilon*(diff(diff(E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t), t), t)))

-E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t)*(k__x^2+k__y^2+k__z^2) = -mu*varepsilon*E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*omega^2*exp(-I*omega*t)

(4.4.8)

-(-E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*exp(-I*omega*t)*(k__x^2+k__y^2+k__z^2) = -mu*varepsilon*E*exp(I*k__x*x)*exp(I*k__y*y)*exp(I*k__z*z)*omega^2*exp(-I*omega*t))/E__wave_(x, y, z, t)

k__x^2+k__y^2+k__z^2 = mu*varepsilon*omega^2

(4.4.9)

So we now have the relationship between k and ω.

 

If we define a wave vector k:

k_ = _i*k__x+_j*k__y+_k*k__z

k_ = _i*k__x+_j*k__y+_k*k__z

(4.4.10)

We can see that:

(k_ = _i*k__x+_j*k__y+_k*k__z).(k_ = _i*k__x+_j*k__y+_k*k__z)

Physics:-Vectors:-Norm(k_)^2 = k__x^2+k__y^2+k__z^2

(4.4.11)

NULL

And so using (4.4.9):

subs(rhs(Physics[Vectors][Norm](k_)^2 = k__x^2+k__y^2+k__z^2) = lhs(Physics[Vectors][Norm](k_)^2 = k__x^2+k__y^2+k__z^2), k__x^2+k__y^2+k__z^2 = mu*varepsilon*omega^2)

Physics:-Vectors:-Norm(k_)^2 = mu*varepsilon*omega^2

(4.4.12)

And so:NULL

LinearAlgebra[Norm](lhs(k_ = _i*k__x+_j*k__y+_k*k__z)) = sqrt(mu*epsilon)*omega

Physics:-Vectors:-Norm(k_) = (mu*varepsilon)^(1/2)*omega

(4.4.13)

Finding the equation for the Electric Field Vector

 

"`B__flux_`(x,y,z,t):='`B__flux_`(x,y,z,t)':`E__field_`(x,y,z,t):='`E__field_`(x,y,z,t)':`H__field_`(x,y,z,t):='`H__field_`(x,y,z,t)':E_:='E_':H_:='H_':k_:='k_':r_:='r_':"

 

Recall Gauss' Law and the first of the Constitutive Relations:

Physics[Vectors][Divergence](D__flux_(x, y, z, t)) = 0

Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)) = 0

(4.5.1)

D__flux_(x, y, z, t) = varepsilon*E__field_(x, y, z, t)

D__flux_(x, y, z, t) = varepsilon*E__field_(x, y, z, t)

(4.5.2)

Substitute D for εE:

subs(D__flux_(x, y, z, t) = varepsilon*E__field_(x, y, z, t), Physics[Vectors][Divergence](D__flux_(x, y, z, t)) = 0)

Physics:-Vectors:-Divergence(varepsilon*E__field_(x, y, z, t)) = 0

(4.5.3)

Construct an expression for the electric field that will satisfy Maxwell's equations:NULL

E__field_(x, y, z, t) = E_*exp(I*(omega*t-k_.r_))

E__field_(x, y, z, t) = E_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_, r_)))

(4.5.4)

Define the wave vector k, position vector r and electric field vector E:``

k_ = _i*k__x+_j*k__y+_k*k__z

k_ = _i*k__x+_j*k__y+_k*k__z

(4.5.5)

r_ = _i*x+_j*y+_k*z

r_ = _i*x+_j*y+_k*z

(4.5.6)

E_ = _i*e__x+_j*e__y+_k*e__z

E_ = _i*e__x+_j*e__y+_k*e__z

(4.5.7)

Substitute in to the expression for electric field:

subs({E_ = _i*e__x+_j*e__y+_k*e__z, k_ = _i*k__x+_j*k__y+_k*k__z, r_ = _i*x+_j*y+_k*z}, E__field_(x, y, z, t) = E_*exp(I*(omega*t-Physics[Vectors][`.`](k_, r_))))

E__field_(x, y, z, t) = (_i*e__x+_j*e__y+_k*e__z)*exp(I*(omega*t-Physics:-Vectors:-`.`(_i*k__x+_j*k__y+_k*k__z, _i*x+_j*y+_k*z)))

(4.5.8)

And substitute into the expression derived from Gauss' Law:

subs(E__field_(x, y, z, t) = (_i*e__x+_j*e__y+_k*e__z)*exp(I*(omega*t-Physics[Vectors][`.`](_i*k__x+_j*k__y+_k*k__z, _i*x+_j*y+_k*z))), Physics[Vectors][Divergence](varepsilon*E__field_(x, y, z, t)) = 0)

varepsilon*Physics:-Vectors:-Divergence((_i*e__x+_j*e__y+_k*e__z)*exp(I*(-k__x*x-k__y*y-k__z*z+omega*t))) = 0

(4.5.9)

Activate the divergence operator:

value(varepsilon*Physics[Vectors][Divergence]((_i*e__x+_j*e__y+_k*e__z)*exp(I*(-k__x*x-k__y*y-k__z*z+omega*t))) = 0)

-I*varepsilon*exp(-I*(k__x*x+k__y*y+k__z*z-omega*t))*(e__x*k__x+e__y*k__y+e__z*k__z) = 0

(4.5.10)

Solve for the components of the E vector:

solve(-I*varepsilon*exp(-I*(k__x*x+k__y*y+k__z*z-omega*t))*(e__x*k__x+e__y*k__y+e__z*k__z) = 0, {e__x, e__y, e__z})

{e__x = -k__y*e__y/k__x-e__z*k__z/k__x, e__y = e__y, e__z = e__z}

(4.5.11)

It appears that the E vector must be perpendicular to the wave vector k.  We can check this as follows:``

k_ := rhs(k_ = _i*k__x+_j*k__y+_k*k__z); r_ := rhs(r_ = _i*x+_j*y+_k*z); E_ := rhs(E_ = _i*e__x+_j*e__y+_k*e__z)

``

solve(E_.k_ = 0, {e__x, e__y, e__z})

{e__x = -(e__y*k__y+e__z*k__z)/k__x, e__y = e__y, e__z = e__z}

(4.5.12)

Yes, since (4.5.11) and (4.5.12) are the same, we have confirmed that the direction of the electric field E, must be perpendicular to the wave vector k - a wave with a transverse electric field.

``

Finding the wave equation for the Magnetic Field Vector

 

 
"`D__flux_`(x,y,z,t):='`D__flux_`(x,y,z,t)':`B__flux_`(x,y,z,t):='`B__flux_`(x,y,z,t)':`E__field_`(x,y,z,t):='`E__field_`(x,y,z,t)':`H__field_`(x,y,z,t):='`H__field_`(x,y,z,t)':E_:='E_':H_:='H_':"

``

We know that the PDEs for the electric and magnetic fields have the same form:``

E__fieldPDE

%Laplacian(E__field_(x, y, z, t)) = mu*varepsilon*(diff(diff(E__field_(x, y, z, t), t), t))

(4.6.1)

H__fieldPDE

%Laplacian(H__field_(x, y, z, t)) = varepsilon*mu*(diff(diff(H__field_(x, y, z, t), t), t))

(4.6.2)

And so we would expect the solution for the magnetic field to have the same form as the solution for the electric field.NULL

 

We know that a solution for  for the electric field equation is:

E__field_ := E_*exp(I*(omega*t-k_.r_))

E_*exp(I*(-k__x*x-k__y*y-k__z*z+omega*t))

(4.6.3)

And so the solution for the magnetic field equation will be:

H__field_ := H_*exp(I*(omega*t-k_.r_))

H_*exp(I*(-k__x*x-k__y*y-k__z*z+omega*t))

(4.6.4)

We used Gauss' Law to prove that the electric field vector must be perpendicular to the wave vector k:

Physics[Vectors][Divergence](D__flux_(x, y, z, t)) = 0

Physics:-Vectors:-Divergence(D__flux_(x, y, z, t)) = 0

(4.6.5)

And we could do the same for the magnetic field.  Recall Gauss' Law for Magnetism and the second of the Constitutive Relations:``

Physics[Vectors][Divergence](B__flux_(x, y, z, t)) = 0

Physics:-Vectors:-Divergence(B__flux_(x, y, z, t)) = 0

(4.6.6)

So the magnetic field vector H must also be perpendicular to the wave vector, k.

````

Proving that E and H are perpendicular

 

"`E__field_`(x, y, z, t):='`E__field_`(x, y, z, t)': `H__field_`(x,y,z,t):='`H__field_`(x,y,z,t)':E_:='E_':H_:='H_':"r_ := 'r_'; k_ := 'k_'``

``

Recall the Maxwell–Faraday equation and the second of the Constitutive Relations:

Maxwell_1 := Curl(E__field_(x, y, z, t)) = -(diff(B__flux_(x, y, z, t), t))

Physics:-Vectors:-Curl(E__field_(x, y, z, t)) = -(diff(B__flux_(x, y, z, t), t))

(4.7.1)

NULL

Eq_2 := B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

(4.7.2)

Substituting H for B:

subs(B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t), Physics[Vectors][Curl](E__field_(x, y, z, t)) = -(diff(B__flux_(x, y, z, t), t)))

Physics:-Vectors:-Curl(E__field_(x, y, z, t)) = -(diff(mu*H__field_(x, y, z, t), t))

(4.7.3)

Construct expressions for the electric field and magnetic field that will satisfy Maxwell's equations:``

E__field_(x, y, z, t) = E_*exp(I*(omega*t-k_(x, y, z).r_(x, y, z)))

E__field_(x, y, z, t) = E_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z))))

(4.7.4)

H__field_(x, y, z, t) = H_*exp(I*(omega*t-k_(x, y, z).r_(x, y, z)))

H__field_(x, y, z, t) = H_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z))))

(4.7.5)

Substitute into (4.7.3):NULL

subs({E__field_(x, y, z, t) = E_*exp(I*(omega*t-Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))), H__field_(x, y, z, t) = H_*exp(I*(omega*t-Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z))))}, Physics[Vectors][Curl](E__field_(x, y, z, t)) = -(diff(mu*H__field_(x, y, z, t), t)))

Physics:-Vectors:-Curl(E_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z))))) = -mu*(diff(H_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))), t))

(4.7.6)

Activate the curl operator:

value(Physics[Vectors][Curl](E_*exp(I*(omega*t-Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z))))) = -mu*(diff(H_*exp(I*(omega*t-Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))), t)))

-I*Physics:-Vectors:-`.`(diff(k_(x, y, z), x), r_(x, y, z))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_i, E_)-I*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), x))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_i, E_)-I*Physics:-Vectors:-`.`(diff(k_(x, y, z), y), r_(x, y, z))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_j, E_)-I*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), y))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_j, E_)-I*Physics:-Vectors:-`.`(diff(k_(x, y, z), z), r_(x, y, z))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_k, E_)-I*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), z))*exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*Physics:-Vectors:-`&x`(_k, E_) = -I*mu*H_*omega*exp(I*(omega*t-Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z))))

(4.7.7)

``

 

And now solve for the H vector:

solve(-I*Physics[Vectors][`.`](diff(k_(x, y, z), x), r_(x, y, z))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_i, E_)-I*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), x))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_i, E_)-I*Physics[Vectors][`.`](diff(k_(x, y, z), y), r_(x, y, z))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_j, E_)-I*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), y))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_j, E_)-I*Physics[Vectors][`.`](diff(k_(x, y, z), z), r_(x, y, z))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_k, E_)-I*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), z))*exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*Physics[Vectors][`&x`](_k, E_) = -I*mu*H_*omega*exp(I*(omega*t-Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))), H_)

exp(I*omega*t-I*Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z)))*(Physics:-Vectors:-`&x`(_i, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), x), r_(x, y, z))+Physics:-Vectors:-`&x`(_i, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), x))+Physics:-Vectors:-`&x`(_j, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), y), r_(x, y, z))+Physics:-Vectors:-`&x`(_j, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), y))+Physics:-Vectors:-`&x`(_k, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), z), r_(x, y, z))+Physics:-Vectors:-`&x`(_k, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), z)))/(exp(-I*(-omega*t+Physics:-Vectors:-`.`(k_(x, y, z), r_(x, y, z))))*mu*omega)

(4.7.8)

simplify(exp(I*omega*t-I*Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z)))*(Physics[Vectors][`&x`](_i, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), x), r_(x, y, z))+Physics[Vectors][`&x`](_i, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), x))+Physics[Vectors][`&x`](_j, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), y), r_(x, y, z))+Physics[Vectors][`&x`](_j, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), y))+Physics[Vectors][`&x`](_k, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), z), r_(x, y, z))+Physics[Vectors][`&x`](_k, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), z)))/(exp(-I*(-omega*t+Physics[Vectors][`.`](k_(x, y, z), r_(x, y, z))))*mu*omega))

(Physics:-Vectors:-`&x`(_i, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), x), r_(x, y, z))+Physics:-Vectors:-`&x`(_i, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), x))+Physics:-Vectors:-`&x`(_j, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), y), r_(x, y, z))+Physics:-Vectors:-`&x`(_j, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), y))+Physics:-Vectors:-`&x`(_k, E_)*Physics:-Vectors:-`.`(diff(k_(x, y, z), z), r_(x, y, z))+Physics:-Vectors:-`&x`(_k, E_)*Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), z)))/(mu*omega)

(4.7.9)

collect((Physics[Vectors][`&x`](_i, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), x), r_(x, y, z))+Physics[Vectors][`&x`](_i, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), x))+Physics[Vectors][`&x`](_j, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), y), r_(x, y, z))+Physics[Vectors][`&x`](_j, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), y))+Physics[Vectors][`&x`](_k, E_)*Physics[Vectors][`.`](diff(k_(x, y, z), z), r_(x, y, z))+Physics[Vectors][`&x`](_k, E_)*Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), z)))/(mu*omega), {`&x`(_i, E_), `&x`(_j, E_), `&x`(_k, E_)})

(Physics:-Vectors:-`.`(diff(k_(x, y, z), x), r_(x, y, z))+Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), x)))*Physics:-Vectors:-`&x`(_i, E_)/(mu*omega)+(Physics:-Vectors:-`.`(diff(k_(x, y, z), y), r_(x, y, z))+Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), y)))*Physics:-Vectors:-`&x`(_j, E_)/(mu*omega)+(Physics:-Vectors:-`.`(diff(k_(x, y, z), z), r_(x, y, z))+Physics:-Vectors:-`.`(k_(x, y, z), diff(r_(x, y, z), z)))*Physics:-Vectors:-`&x`(_k, E_)/(mu*omega)

(4.7.10)

Substitute for the position vector r_ and the wave vector k_:

subs(r_(x, y, z) = _i*x+_j*y+_k*z, k_(x, y, z) = _i*k__x+_j*k__y+_k*k__z, (Physics[Vectors][`.`](diff(k_(x, y, z), x), r_(x, y, z))+Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), x)))*Physics[Vectors][`&x`](_i, E_)/(mu*omega)+(Physics[Vectors][`.`](diff(k_(x, y, z), y), r_(x, y, z))+Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), y)))*Physics[Vectors][`&x`](_j, E_)/(mu*omega)+(Physics[Vectors][`.`](diff(k_(x, y, z), z), r_(x, y, z))+Physics[Vectors][`.`](k_(x, y, z), diff(r_(x, y, z), z)))*Physics[Vectors][`&x`](_k, E_)/(mu*omega))

(Physics:-Vectors:-`.`(diff(_i*k__x+_j*k__y+_k*k__z, x), _i*x+_j*y+_k*z)+Physics:-Vectors:-`.`(_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, x)))*Physics:-Vectors:-`&x`(_i, E_)/(mu*omega)+(Physics:-Vectors:-`.`(diff(_i*k__x+_j*k__y+_k*k__z, y), _i*x+_j*y+_k*z)+Physics:-Vectors:-`.`(_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, y)))*Physics:-Vectors:-`&x`(_j, E_)/(mu*omega)+(Physics:-Vectors:-`.`(diff(_i*k__x+_j*k__y+_k*k__z, z), _i*x+_j*y+_k*z)+Physics:-Vectors:-`.`(_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, z)))*Physics:-Vectors:-`&x`(_k, E_)/(mu*omega)

(4.7.11)

H_ := simplify((Physics[Vectors][`.`](diff(_i*k__x+_j*k__y+_k*k__z, x), _i*x+_j*y+_k*z)+Physics[Vectors][`.`](_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, x)))*Physics[Vectors][`&x`](_i, E_)/(mu*omega)+(Physics[Vectors][`.`](diff(_i*k__x+_j*k__y+_k*k__z, y), _i*x+_j*y+_k*z)+Physics[Vectors][`.`](_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, y)))*Physics[Vectors][`&x`](_j, E_)/(mu*omega)+(Physics[Vectors][`.`](diff(_i*k__x+_j*k__y+_k*k__z, z), _i*x+_j*y+_k*z)+Physics[Vectors][`.`](_i*k__x+_j*k__y+_k*k__z, diff(_i*x+_j*y+_k*z, z)))*Physics[Vectors][`&x`](_k, E_)/(mu*omega))

(k__x*Physics:-Vectors:-`&x`(_i, E_)+k__y*Physics:-Vectors:-`&x`(_j, E_)+k__z*Physics:-Vectors:-`&x`(_k, E_))/(mu*omega)

(4.7.12)

This looks a lot like k cross E.NULL

`&x`(k_, E_)/(mu*omega)

Physics:-Vectors:-`&x`(k_, E_)/(mu*omega)

(4.7.13)

Lets test this by substituting for the wave vector k_:

subs(k_ = _i*k__x+_j*k__y+_k*k__z, Physics[Vectors][`&x`](k_, E_)/(mu*omega))

Physics:-Vectors:-`&x`(_i*k__x+_j*k__y+_k*k__z, E_)/(mu*omega)

(4.7.14)

expand(Physics[Vectors][`&x`](_i*k__x+_j*k__y+_k*k__z, E_)/(mu*omega))

k__x*Physics:-Vectors:-`&x`(_i, E_)/(mu*omega)+k__y*Physics:-Vectors:-`&x`(_j, E_)/(mu*omega)+k__z*Physics:-Vectors:-`&x`(_k, E_)/(mu*omega)

(4.7.15)

possibleH_ := collect(k__x*Physics[Vectors][`&x`](_i, E_)/(mu*omega)+k__y*Physics[Vectors][`&x`](_j, E_)/(mu*omega)+k__z*Physics[Vectors][`&x`](_k, E_)/(mu*omega), {mu, omega})

(k__x*Physics:-Vectors:-`&x`(_i, E_)+k__y*Physics:-Vectors:-`&x`(_j, E_)+k__z*Physics:-Vectors:-`&x`(_k, E_))/(mu*omega)

(4.7.16)

We were correct:

H_-possibleH_

0

(4.7.17)

So we have shown that:

H_ := Physics[Vectors][`&x`](k_, E_)/(mu*omega)

Physics:-Vectors:-`&x`(k_, E_)/(mu*omega)

(4.7.18)

And so:

H__field_(x, y, z, t) = `&x`(k_, E__field_(x, y, z, t))/(mu*omega)

H__field_(x, y, z, t) = Physics:-Vectors:-`&x`(k_, E__field_(x, y, z, t))/(mu*omega)

(4.7.19)

Since we already know from (4.5.11) that the electric field vector E_ is perpendicular to the wave vector k_, then (4.7.18) proves that E_ and H_ must also be perpendicular.NULL

``NULL

Calculating wave impedance

 

E_ := 'E_'; H_ := 'H_'; k_ := 'k_'Z := 'Z'

``

We know from (4.7.18) that the magnetic field vector H_ and the electric field vector E_ are related as follows:

H_ = Physics[Vectors][`&x`](k_, E_)/(mu*omega)

H_ = Physics:-Vectors:-`&x`(k_, E_)/(mu*omega)

(5.1)

And so:

LinearAlgebra[Norm](H_) = LinearAlgebra[Norm](k_)*LinearAlgebra[Norm](E_)/(mu*omega)

Physics:-Vectors:-Norm(H_) = Physics:-Vectors:-Norm(k_)*Physics:-Vectors:-Norm(E_)/(mu*omega)

(5.2)

The wave impedance Z is defined as the ratio between the magnitude of the electric field ||E_|| and the magnitude of the magnetic field ||H_||:

Z = lhs(((Physics[Vectors][Norm](H_) = Physics[Vectors][Norm](k_)*Physics[Vectors][Norm](E_)/(mu*omega))*mu)*omega/(LinearAlgebra[Norm](H_)*LinearAlgebra[Norm](k_)))

Z = mu*omega/Physics:-Vectors:-Norm(k_)

(5.3)

We know from (4.4.13) that:

Physics[Vectors][Norm](k_) = (mu*varepsilon)^(1/2)*omega

Physics:-Vectors:-Norm(k_) = (mu*varepsilon)^(1/2)*omega

(5.4)

So we can substitute for ||k_||:

subs(Physics[Vectors][Norm](k_) = (mu*varepsilon)^(1/2)*omega, Z = mu*omega/Physics[Vectors][Norm](k_))

Z = mu/(mu*varepsilon)^(1/2)

(5.5)

And so:

Z = sqrt(mu/epsilon)

Z = (mu/varepsilon)^(1/2)

(5.6)

(5.6)If we assume that we are dealing with waves propagating through free space:NULL

`ϵ__0` := evalf(ScientificConstants:-Constant(epsilon[0], units))

0.8854187815e-11*Units:-Unit(A^2*s^4/(m^3*kg))

(5.7)

`μ__0` := evalf(ScientificConstants:-Constant(mu[0], units))

0.1256637062e-5*Units:-Unit(m*kg/(A^2*s^2))

(5.8)

And so the wave impedance will be:

simplify(sqrt(`μ__0`/`ϵ__0`))

376.7303136*Units:-Unit(`Ω`)

(5.9)

As an engineer, I have always thought that the wave impedance in free space was 120π, but it seems that this is not quite true:

evalf(120*Pi)

376.9911185

(5.10)

Calculating wavelength

 

r_ := 'r_'; k_ := 'k_'; E_ := 'E_'NULL``

``

We know that the general solution for the monochromatic (single-frequency) electric field is:``

E__field_ := E_*exp(I*(omega*t-k_.r_))

E_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_, r_)))

(6.1)

The wavelength is defined as the distance between two points in space that have the same field strength.

subs(r_ = r__2*_, t = t__2, E_*exp(I*(omega*t-Physics[Vectors][`.`](k_, r_)))) = subs(r_ = r__1*_, t = t__1, E_*exp(I*(omega*t-Physics[Vectors][`.`](k_, r_))))

E_*exp(I*(omega*t__2-Physics:-Vectors:-`.`(k_, r__2*_))) = E_*exp(I*(omega*t__1-Physics:-Vectors:-`.`(k_, r__1*_)))

(6.2)

And so:

omega*t__2-k_.r__2_ = omega*t__1-k_.r__1_

omega*t__2-Physics:-Vectors:-`.`(k_, r__2_) = omega*t__1-Physics:-Vectors:-`.`(k_, r__1_)

(6.3)

Which can be re-arranged as follows:

omega*(t__2-t__1) = k_.(r__2_-r__1_)

omega*(t__2-t__1) = Physics:-Vectors:-`.`(k_, r__2_-r__1_)

(6.4)

Since we know that the two points will be 2π radians apart:

subs(omega*(t__2-t__1) = 2*Pi, omega*(t__2-t__1) = Physics[Vectors][`.`](k_, r__2_-r__1_))

2*Pi = Physics:-Vectors:-`.`(k_, r__2_-r__1_)

(6.5)

And since we have defined the wavelength as the distance between two points in space that have the same field strength:

subs(k_.(r__2_-r__1_) = LinearAlgebra[Norm](k_)*lambda, 2*Pi = Physics[Vectors][`.`](k_, r__2_-r__1_))

2*Pi = Physics:-Vectors:-Norm(k_)*lambda

(6.6)

And so:

WaveLength := lambda = solve(2*Pi = Physics[Vectors][Norm](k_)*lambda, lambda)

lambda = 2*Pi/Physics:-Vectors:-Norm(k_)

(6.7)

``

Calculating wave velocity

 

"r_(t):='r_(t)': k_:='k_': E_:='E_':"

``

We know that the general solution for the monochromatic (single-frequency) electric field is:NULL

E__field_ := E_*exp(I*(omega*t-k_.r_(t)))

E_*exp(I*(omega*t-Physics:-Vectors:-`.`(k_, r_(t))))

(7.1)

The wave velocity is the speed at which the increase in phase due to the ω*t term is cancelled out by a corresponding increase in the k.r term.  That is:

omega*t-k_.r_(t) = constant

omega*t-Physics:-Vectors:-`.`(k_, r_(t)) = constant

(7.2)

Differentiating both sides with respect to time:

diff(omega*t-Physics[Vectors][`.`](k_, r_(t)) = constant, t)

omega-Physics:-Vectors:-`.`(k_, diff(r_(t), t)) = 0

(7.3)

And so:

v = omega/LinearAlgebra[Norm](k_)

v = omega/Physics:-Vectors:-Norm(k_)

(7.4)

Since we know that:

Physics[Vectors][Norm](k_) = (mu*varepsilon)^(1/2)*omega

Physics:-Vectors:-Norm(k_) = (mu*varepsilon)^(1/2)*omega

(7.5)

WaveVelocity := subs(Physics[Vectors][Norm](k_) = (mu*varepsilon)^(1/2)*omega, v = omega/Physics[Vectors][Norm](k_))

v = 1/(mu*varepsilon)^(1/2)

(7.6)

If we assume that we are dealing with waves propagating through free space:``

`ϵ__0` := evalf(ScientificConstants:-Constant(epsilon[0], units))

0.8854187815e-11*Units:-Unit(A^2*s^4/(m^3*kg))

(7.7)

`μ__0` := evalf(ScientificConstants:-Constant(mu[0], units))

0.1256637062e-5*Units:-Unit(m*kg/(A^2*s^2))

(7.8)

And so the wave velocity will be:

simplify(1/sqrt(`μ__0`*`ϵ__0`))

299792458.0*Units:-Unit(m/s)

(7.9)

Which is the speed of light:

evalf(ScientificConstants:-Constant(c, units))

299792458.*Units:-Unit(m/s)

(7.10)

299792458.0*Units:-Unit(m/s)/(299792458.*Units:-Unit(m/s))

1.000000000

(7.11)

So we have shown that:

c = 1/sqrt(mu*epsilon)

c = 1/(mu*varepsilon)^(1/2)

(7.12)

 

Calculating the ratio between the electric field strength and the magnetic flux density

 

We know that the ratio between the magnitude of the electric field strength vector E_ and the magnetic field strength vector H_ is given by:

E/H = (Z = (mu/varepsilon)^(1/2))

E/H = (Z = (mu/varepsilon)^(1/2))

(8.1)

And we know that the magnetic field strength H_ is related to magnetic flux density as follows:

B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

B__flux_(x, y, z, t) = mu*H__field_(x, y, z, t)

(8.2)

So substituting:

E*mu/B = rhs(E/H = (Z = (mu/varepsilon)^(1/2)))

E*mu/B = (Z = (mu/varepsilon)^(1/2))

(8.3)

E/B = 1/sqrt(mu*epsilon)

E/B = 1/(mu*varepsilon)^(1/2)

(8.4)

And so:

subs(rhs(c = 1/(mu*varepsilon)^(1/2)) = lhs(c = 1/(mu*varepsilon)^(1/2)), E/B = 1/(mu*varepsilon)^(1/2))

E/B = c

(8.5)

Amusing how this boils down to such a simple result :-)

 

The Poynting Vector

 

Defined as S=ExH

 

"`E__field_`(x, y, z, t):='`E__field_`(x, y, z, t)': `H__field_`(x,y,z,t):='`H__field_`(x,y,z,t)':"

``

The Poynting vector represents the directional energy flux (the energy transfer per unit area per unit time) of an electromagnetic field and is defined as follows:

S_(x, y, z, t) = `&x`(E__field_(x, y, z, t), H__field_(x, y, z, t))

S_(x, y, z, t) = Physics:-Vectors:-`&x`(E__field_(x, y, z, t), H__field_(x, y, z, t))

(9.1.1)

This is analogous to a resistive circuit where power is the product of voltage and current: P=V*I.

 

We derived an expression for the magnetic field expressed in terms of the electric field:

H__field_(x, y, z, t) = Physics[Vectors][`&x`](k_, E__field_(x, y, z, t))/(mu*omega)

H__field_(x, y, z, t) = Physics:-Vectors:-`&x`(k_, E__field_(x, y, z, t))/(mu*omega)

(9.1.2)

And we can use this in order to express the Poynting Vector in terms of the electric field only:

subs(H__field_(x, y, z, t) = Physics[Vectors][`&x`](k_, E__field_(x, y, z, t))/(mu*omega), S_(x, y, z, t) = Physics[Vectors][`&x`](E__field_(x, y, z, t), H__field_(x, y, z, t)))

S_(x, y, z, t) = Physics:-Vectors:-`&x`(E__field_(x, y, z, t), Physics:-Vectors:-`&x`(k_, E__field_(x, y, z, t))/(mu*omega))

(9.1.3)

``

Using the vector identity Ax(BxC)=(A.C)B-(A.B)C, this becomes:

S_(x, y, z, t) = ((E__field_(x, y, z, t).E__field_(x, y, z, t))*k_-(E__field_(x, y, z, t).k_)*E__field_(x, y, z, t))/(mu*omega)

S_(x, y, z, t) = (Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*k_-Physics:-Vectors:-`.`(E__field_(x, y, z, t), k_)*E__field_(x, y, z, t))/(mu*omega)

(9.1.4)

Since we know from (4.5.11) that the electric field vector is perpendicular to the wave vector, the second term must be zero, and so:

S_(x, y, z, t) = (E__field_(x, y, z, t).E__field_(x, y, z, t))*k_/(mu*omega)

S_(x, y, z, t) = Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*k_/(mu*omega)

(9.1.5)

We can split the wave vector into a scalar, magnitude component ||k|| and a directional, unit vector component k/||k|| as follows k=||k||*k/||k||, and since we know that:

Physics[Vectors][Norm](k_) = (mu*varepsilon)^(1/2)*omega

Physics:-Vectors:-Norm(k_) = (mu*varepsilon)^(1/2)*omega

(9.1.6)

The :Poynting Vector becomes:

subs(k_ = sqrt(mu*epsilon)*omega*k_/LinearAlgebra[Norm](k_), S_(x, y, z, t) = Physics[Vectors][Norm](E__field_(x, y, z, t))^2*k_/(mu*omega))

S_(x, y, z, t) = Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*(mu*varepsilon)^(1/2)*k_/(Physics:-Vectors:-Norm(k_)*mu)

(9.1.7)

S_(x, y, z, t) = sqrt(epsilon/mu)*(E__field_(x, y, z, t).E__field_(x, y, z, t))*k_/LinearAlgebra[Norm](k_)

S_(x, y, z, t) = (varepsilon/mu)^(1/2)*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*k_/Physics:-Vectors:-Norm(k_)

(9.1.8)

We have already defined wave impedance, Z in terms of μ and ε:``

Z = (mu/varepsilon)^(1/2)

Z = (mu/varepsilon)^(1/2)

(9.1.9)

S_(x, y, z, t) = (E__field_(x, y, z, t).E__field_(x, y, z, t))*k_/(Z*LinearAlgebra[Norm](k_))

S_(x, y, z, t) = Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*k_/(Z*Physics:-Vectors:-Norm(k_))

(9.1.10)

So we have shown that the Poynting Vector points in the direction of the wave vector k and has an instantaneous value equal to the square of the magnitude of the electric field divided by the wave impedance.  This is analogous to a resistive circuit where power, P=V2/R.

 

Since the The wave impedance Z is defined as the ratio between the magnitude of the electric field and the magnitude of the magnetic field, we can also express the Poynting Vector in terms of the magnetic field:

subs(E__field_(x, y, z, t) = Z*H__field_(x, y, z, t), S_(x, y, z, t) = Physics[Vectors][Norm](E__field_(x, y, z, t))^2*k_/(Z*Physics[Vectors][Norm](k_)))

S_(x, y, z, t) = Physics:-Vectors:-Norm(Z*H__field_(x, y, z, t))^2*k_/(Z*Physics:-Vectors:-Norm(k_))

(9.1.11)

S_(x, y, z, t) = Z*(H__field_(x, y, z, t).H__field_(x, y, z, t))*k_/LinearAlgebra[Norm](k_)

S_(x, y, z, t) = Z*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2*k_/Physics:-Vectors:-Norm(k_)

(9.1.12)

This is analogous to a resistive circuit where power, P=I2R.````

Derivation based on the E and H field energy densities

 

"`E__field_`(x, y, z, t):='`E__field_`(x, y, z, t)': `H__field_`(x,y,z,t):='`H__field_`(x,y,z,t)':"

NULL

If we accept that the electric field energy density is given by:

w__e = (1/2)*epsilon*LinearAlgebra[Norm](E__field_(x, y, z, t))^2

w__e = (1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2

(9.2.1)

And that the magnetic field energy density is given by:

w__h = (1/2)*mu*LinearAlgebra[Norm](H__field_(x, y, z, t))^2

w__h = (1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2

(9.2.2)

Then the energy density of the electromagnetic field will be:

w = w__e+w__h

w = w__e+w__h

(9.2.3)

subs(w__e = (1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2, w__h = (1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2, w = w__e+w__h)

w = (1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2

(9.2.4)

The energy contained in a small volume, dV, will be:

`δu` = rhs(w = (1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*`δV`

`δu` = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*`δV`

(9.2.5)

If we define the volume as having a surface with area dA and a length of dl:

subs(`δV` = `δA`*`δl`, `δu` = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*`δV`)

`δu` = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*`δA`*`δl`

(9.2.6)

We now assume that the volume is oriented in space such that the surface dA is normal to the wave vector k and the length dl is parallel to k.  Since we know that the electromagnetic wave propagates in the direction of k, this means that all of the energy contained in the volume dV will exit through the surface dA at the speed of light, c.  The time required for this to happen will be dt=dl/c seconds.  We can therefore replace dl with c*dt:

subs(`δl` = c*`δt`, `δu` = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*`δA`*`δl`)

`δu` = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*`δA`*c*`δt`

(9.2.7)

Power is defined as the rate of doing work, so we can calculate power by dividing du by dt:

``

(`δu` = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*`δA`*c*`δt`)*(1/`δt`)

`δu`/`δt` = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*`δA`*c

(9.2.8)

And power flux density is defied as power per unit area:

(`δu`/`δt` = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*`δA`*c)*(1/`δA`)

`δu`/(`δt`*`δA`) = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*c

(9.2.9)

Since we know that power is flowing in the direction of the wave vector k, we can express this relation as a vector quantity:

S_(x, y, z, t) = rhs(`δu`/(`δA`*`δt`) = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*c)*k_/LinearAlgebra[Norm](k_)

S_(x, y, z, t) = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*c*k_/Physics:-Vectors:-Norm(k_)

(9.2.10)

This should be the Poynting Vector.

 

We know that:

c = 1/(mu*varepsilon)^(1/2)

c = 1/(mu*varepsilon)^(1/2)

(9.2.11)

And so:

subs(c = 1/(mu*varepsilon)^(1/2), S_(x, y, z, t) = ((1/2)*varepsilon*Physics[Vectors][Norm](E__field_(x, y, z, t))^2+(1/2)*mu*Physics[Vectors][Norm](H__field_(x, y, z, t))^2)*c*k_/Physics[Vectors][Norm](k_))

S_(x, y, z, t) = ((1/2)*varepsilon*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(1/2)*mu*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*k_/((mu*varepsilon)^(1/2)*Physics:-Vectors:-Norm(k_))

(9.2.12)

Which can be re-arranged:

S_(x, y, z, t) = (1/2*(sqrt(epsilon/mu)*LinearAlgebra[Norm](E__field_(x, y, z, t))^2+sqrt(mu/epsilon)*LinearAlgebra[Norm](H__field_(x, y, z, t))^2))*k_/LinearAlgebra[Norm](k_)

S_(x, y, z, t) = (1/2)*((varepsilon/mu)^(1/2)*Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2+(mu/varepsilon)^(1/2)*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*k_/Physics:-Vectors:-Norm(k_)

(9.2.13)

We know that:

Z = (mu/varepsilon)^(1/2)

Z = (mu/varepsilon)^(1/2)

(9.2.14)

And so:

S_(x, y, z, t) = (1/2*(LinearAlgebra[Norm](E__field_(x, y, z, t))^2/Z+Z*LinearAlgebra[Norm](H__field_(x, y, z, t))^2))*k_/LinearAlgebra[Norm](k_)

S_(x, y, z, t) = (1/2)*(Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2/Z+Z*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2)*k_/Physics:-Vectors:-Norm(k_)

(9.2.15)

If we compare this result with:

S_(x, y, z, t) = Physics[Vectors][Norm](E__field_(x, y, z, t))^2*k_/(Z*Physics[Vectors][Norm](k_))

S_(x, y, z, t) = Physics:-Vectors:-Norm(E__field_(x, y, z, t))^2*k_/(Z*Physics:-Vectors:-Norm(k_))

(9.2.16)

And:

S_(x, y, z, t) = Z*Physics[Vectors][Norm](H__field_(x, y, z, t))^2*k_/Physics[Vectors][Norm](k_)

S_(x, y, z, t) = Z*Physics:-Vectors:-Norm(H__field_(x, y, z, t))^2*k_/Physics:-Vectors:-Norm(k_)

(9.2.17)

We see that all three expressions for S are equivalent, and can also see that the power is split evenly between the electric field and magnetic field components of the electromagnetic wave :-)

````

Conservation of energy

 

As the electromagnetic wave moves forward, the surface area of the wave will increase as the square of the distance traveled.  Since we are considering the propagation of electromagnetic waves in a vacuum, the energy in the wave will remain constant with no energy being lost to heating effects.  This means that the Poynting Vector must reduce as the square of the distance traveled.  Using (9.2.16) and (9.2.17) above, we see that the magnitude of the electric and magnetic fields must therefore reduce linearly with distance.``

``

Visualising the electric field

 

Procedures

   

Plane waves

 

The electric field strength is calculated for a plane wave (with a constant wave vector) in the XY plane at height z.

Explore(PlaneWaveEfieldDensityPlotInXYplane(z, k__x, k__y, k__z, omega, t), parameters = [[z = -5 .. 5, animate = false], [k__x = -5 .. 5, animate = false], [k__y = -5 .. 5, animate = false], [k__z = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [z = 1, k__x = 1, k__y = 1, k__z = 1, omega = .1, t = 0], title = "Electric Field", animate = true)

The electric field strength is calculated for a plane wave (with a constant wave vector) in the XY plane at height z.

Explore(PlaneWaveEfieldPlot3dInXYplane(z, k__x, k__y, k__z, omega, t), parameters = [[z = -5 .. 5, animate = false], [k__x = -5 .. 5, animate = false], [k__y = -5 .. 5, animate = false], [k__z = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [z = 1, k__x = 1, k__y = 1, k__z = 1, omega = .1, t = 0], title = "Electric Field", animate = true)

The electric field vector is calculated for a plane wave (with a constant wave vector).  It is easier to examine the electric field by rotating the plot without running the animation.

Explore(PlaneWaveEfieldPlot(k__x, k__y, k__z, omega, t), parameters = [[k__x = -5 .. 5, animate = false], [k__y = -5 .. 5, animate = false], [k__z = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [k__x = 2, k__y = 0, k__z = 0, omega = .1, t = 0], title = "Electric Field", animate = true)

NULL

Spherical waves

 

The electric field strength is calculated for a spherical wave (with the wave vector having a fixed magnitude and pointing out radially from the origin) in the XY plane at height z.  In order for energy to be conserved, the amplitude of the wave should reduce linearly with distance, but this has not been done in the interest of better aesthetics.

Explore(SphericalWaveEfieldDensityPlotInXYplane(z, k, omega, t), parameters = [[z = -10 .. 10, animate = false], [k = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [z = 1, k = 2, omega = .1, t = 0], title = "Electric Field", animate = true)

The electric field strength is calculated for a spherical wave (with the wave vector having a fixed magnitude and pointing out radially from the origin) in the XY plane at height z.  It is interesting to watch the spherical wave become closer to a plane wave by increasing z.  In order for energy to be conserved, the amplitude of the wave should reduce linearly with distance, but this has not been done in the interest of better aesthetics.

Explore(SphericalWaveEfieldPlot3dInXYplane(z, k, omega, t), parameters = [[z = -10 .. 10, animate = false], [k = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [z = 1, k = 2, omega = .1, t = 0], title = "Electric Field", animate = true)

The electric field vector is calculated for a spherical wave (with the wave vector having a fixed magnitude and pointing out radially from the origin).  It is easier to examine the electric field by rotating the plot without running the animation.  In order for energy to be conserved, the amplitude of the wave should reduce linearly with distance, but this has not been done in the interest of better aesthetics.

Explore(SphericalWaveEfieldPlot(k, omega, t), parameters = [[k = -5 .. 5, animate = false], [omega = .1 .. 1, animate = false], [t = 0 .. 100, animate = true]], initialvalues = [k = 2, omega = .1, t = 0], title = "Electric Field", animate = true)

NULLNULL

NULL````


 

Download Examining_Maxwells_Equations.mw

@Pascal4QM Thank you Pascal, your setting look a lot better :-)

I have managed to show why E and H must be perpendicular, so here is the new vesion...

Calculating the wave impedance is still work in progress.

Steve

Examining_Maxwells_Equations.mw

@acer Thank you acer, this all makes perfect sense now.

The unevaluation quotes make a really useful tool :-). And I had not appreciated the subtlety of how they work when used to unassign variables.

Also thanks for "if not [x,y]::list(numeric) then return 'procname'(args); end if;", I would have never thought of doing that!

 

 

@vv Thank you for this.

As you say, this appears to be a more general solution because it even works with my clumsy PUV procedure.  But I can't work out how it works?

 

@Preben Alsholm Thank you for introducing me to piecewise, it works well with fieldplot3d!

@Kitonum Thank you Kitonum,

Your versions of UV and PUV are really nice :-)

I didn't realise that arctan would accept two-arguments, but this is a really nice way of avoiding my clumsy if statements!

Thanks for helping.

 

@Carl Love Thank you Carl!

Can I ask where to find help on this please?: x:= `tools/genglobal`(:-_x)

Also, what is the correct format for the second argument is this version of PV?

 

@Kitonum Thank you Kitonum,

As a Maple newbie, I had not appreciated the possibility of solving for more than one variable ({x,y,z} vs {x})!

By the way, do you know a way to plot an arrow in order to visualise the vector <a,b,c> protruding from the perpendicular vector plane?

Finding a perpendicular vector - general solution

restart

with(plots)

NULL

 Kitonum's general solution (as a Maple newbie, I had not appreciated the possibility of solving for more than one variable ({x,y,z} vs {x}):

PV := proc (a, b, c) options operator, arrow; solve(a*x+b*y+c*z = 0, {x, y, z}) end proc

proc (a, b, c) options operator, arrow; solve(a*x+b*y+c*z = 0, {x, y, z}) end proc

(1)

PV(a, b, c)

{x = -(b*y+c*z)/a, y = y, z = z}

(2)

PV(0, b, c)

{x = x, y = -c*z/b, z = z}

(3)

PV(a, 0, c)

{x = -c*z/a, y = y, z = z}

(4)

PV(a, b, 0)

{x = -b*y/a, y = y, z = z}

(5)

PV(0, 0, c)

{x = x, y = y, z = 0}

(6)

PV(0, b, 0)

{x = x, y = 0, z = z}

(7)

PV(a, 0, 0)

{x = 0, y = y, z = z}

(8)

PV(1, 2, 3)

{x = -2*y-3*z, y = y, z = z}

(9)

So it always seems to work :-)

We can visualise the perpendicular vector plane like this:

Explore(implicitplot3d(a*x+b*y+c*z = 0, x = -3 .. 3, y = -3 .. 3, z = -3 .. 3), parameters = [a = -2 .. 2, b = -2 .. 2, c = -2 .. 2], initialvalues = [a = 1, b = 0, c = 0])

NULL

 

Download Finding_a_perpendicular_vector_-_general_solution.mw

@ecterrab Thank you Edgardo,

I attempted another method using Rodrigues' rotation formula, with partial success:

 

Finding a Perpendicular Vector using Rodrigues' rotation formula

restart

with(Physics)with(Vectors)

[`&x`, `+`, `.`, ChangeBasis, ChangeCoordinates, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, Setup, diff]

(1)

Setup(mathematicalnotation = true)

[mathematicalnotation = true]

(2)

Define a vectors k:

k_ := _i*k__x+_j*k__y+_k*k__z

_i*k__x+_j*k__y+_k*k__z

(3)

We wish to find a vector p that is perpendicular to k (p for perpendicular) and can achieve this using the following three step method:

1. 

Rotate the vector k around the x-axis until the z component is zero.

2. 

Create a vector that is perpendicular to the rotated vector by setting x_new = -y_old and y_new = x_old.

3. 

Rotate the vector back again, to create a vector p that is perpendicular to k.

 

We use Rodrigues' rotation formula:

Rodrigues := Vrot_ = V_*cos(angle)+`&x`(axis_, V_)*sin(angle)+axis_*(axis_.V_)*(1-cos(angle))

Vrot_ = V_*cos(angle)+Physics:-Vectors:-`&x`(axis_, V_)*sin(angle)-axis_*Physics:-Vectors:-`.`(axis_, V_)*cos(angle)+axis_*Physics:-Vectors:-`.`(axis_, V_)

(4)

Now we find the rotation angle that will deliver a _k component of zero.  We will rotate about the _i axis:

subs(V_ = k_, axis_ = _i, Vrot_ = V_*cos(angle)+Physics[Vectors][`&x`](axis_, V_)*sin(angle)-axis_*Physics[Vectors][`.`](axis_, V_)*cos(angle)+axis_*Physics[Vectors][`.`](axis_, V_))

Vrot_ = (_i*k__x+_j*k__y+_k*k__z)*cos(angle)+Physics:-Vectors:-`&x`(_i, _i*k__x+_j*k__y+_k*k__z)*sin(angle)-_i*Physics:-Vectors:-`.`(_i, _i*k__x+_j*k__y+_k*k__z)*cos(angle)+_i*Physics:-Vectors:-`.`(_i, _i*k__x+_j*k__y+_k*k__z)

(5)

For the _k component to be zero, the dot product of the rotated vector Vrot_ and _k must be zero:

rhs(Vrot_ = (_i*k__x+_j*k__y+_k*k__z)*cos(angle)+Physics[Vectors][`&x`](_i, _i*k__x+_j*k__y+_k*k__z)*sin(angle)-_i*Physics[Vectors][`.`](_i, _i*k__x+_j*k__y+_k*k__z)*cos(angle)+_i*Physics[Vectors][`.`](_i, _i*k__x+_j*k__y+_k*k__z))._k = 0

cos(angle)*k__z+sin(angle)*k__y = 0

(6)

solve({cos(angle)*k__z+sin(angle)*k__y = 0}, angle)

{angle = -arctan(k__z/k__y)}

(7)

NULL

OK, now we have the angle, we can go ahead and rotate vector k:

subs(V_ = k_, axis_ = _i, {angle = -arctan(k__z/k__y)}, Vrot_ = V_*cos(angle)+Physics[Vectors][`&x`](axis_, V_)*sin(angle)-axis_*Physics[Vectors][`.`](axis_, V_)*cos(angle)+axis_*Physics[Vectors][`.`](axis_, V_))

Vrot_ = (_i*k__x+_j*k__y+_k*k__z)*cos(-arctan(k__z/k__y))+Physics:-Vectors:-`&x`(_i, _i*k__x+_j*k__y+_k*k__z)*sin(-arctan(k__z/k__y))-_i*Physics:-Vectors:-`.`(_i, _i*k__x+_j*k__y+_k*k__z)*cos(-arctan(k__z/k__y))+_i*Physics:-Vectors:-`.`(_i, _i*k__x+_j*k__y+_k*k__z)

(8)

simplify(Vrot_ = (_i*k__x+_j*k__y+_k*k__z)*cos(-arctan(k__z/k__y))+Physics[Vectors][`&x`](_i, _i*k__x+_j*k__y+_k*k__z)*sin(-arctan(k__z/k__y))-_i*Physics[Vectors][`.`](_i, _i*k__x+_j*k__y+_k*k__z)*cos(-arctan(k__z/k__y))+_i*Physics[Vectors][`.`](_i, _i*k__x+_j*k__y+_k*k__z))

Vrot_ = (k__x*_i*((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y+_j*(k__y^2+k__z^2))/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y)

(9)

Now we create a vector Prot which is perpendicular to Vrot.  Since there is no _k component, we have reduced this to a simple 2D problem:

oldX := rhs(Vrot_ = (k__x*_i*((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y+_j*(k__y^2+k__z^2))/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y))._i

k__x

(10)

oldY := rhs(Vrot_ = (k__x*_i*((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y+_j*(k__y^2+k__z^2))/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y))._j

k__y/(k__z^2/k__y^2+1)^(1/2)+k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y)

(11)

Prot_ := -_i*oldY+_j*oldX

_i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j

(12)

Check that Vrot and Prot are perpendicular:

rhs(Vrot_ = (k__x*_i*((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y+_j*(k__y^2+k__z^2))/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y)).Prot_

0

(13)

NULL

Now we rotate the Prot vector back again to produce a vector p that should be perpendicular to k:

subs(V_ = Prot_, axis_ = _i, angle = arctan(k__z/k__y), Vrot_ = V_*cos(angle)+Physics[Vectors][`&x`](axis_, V_)*sin(angle)-axis_*Physics[Vectors][`.`](axis_, V_)*cos(angle)+axis_*Physics[Vectors][`.`](axis_, V_))

Vrot_ = (_i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*cos(arctan(k__z/k__y))+Physics:-Vectors:-`&x`(_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*sin(arctan(k__z/k__y))-_i*Physics:-Vectors:-`.`(_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*cos(arctan(k__z/k__y))+_i*Physics:-Vectors:-`.`(_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)

(14)

simplify(Vrot_ = (_i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*cos(arctan(k__z/k__y))+Physics[Vectors][`&x`](_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*sin(arctan(k__z/k__y))-_i*Physics[Vectors][`.`](_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j)*cos(arctan(k__z/k__y))+_i*Physics[Vectors][`.`](_i, _i*(-k__y/(k__z^2/k__y^2+1)^(1/2)-k__z^2/((k__z^2/k__y^2+1)^(1/2)*k__y))+k__x*_j))

Vrot_ = -(_i*k__y^2+_i*k__z^2-_j*k__x*k__y-_k*k__x*k__z)/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y)

(15)

p_ := rhs(Vrot_ = -(_i*k__y^2+_i*k__z^2-_j*k__x*k__y-_k*k__x*k__z)/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y))

-(_i*k__y^2+_i*k__z^2-_j*k__x*k__y-_k*k__x*k__z)/(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y)

(16)

Let's see is p is indeed perpendicular to k:

k_.p_

0

(17)

It works!

 

It looks as though we should be able to re-arrange p like this (assigned to q to preserve the original p):

q_ := -(k__y^2+k__z^2)*_i/sqrt(k__y^2+k__z^2)+k__x*k__y*_j/sqrt(k__y^2+k__z^2)+k__x*k__z*_k/sqrt(k__y^2+k__z^2)

-(k__y^2+k__z^2)^(1/2)*_i+k__x*k__y*_j/(k__y^2+k__z^2)^(1/2)+k__x*k__z*_k/(k__y^2+k__z^2)^(1/2)

(18)

But something is wrong, because q is not the same as p:

simplify(q_-p_)

-(((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y-(k__y^2+k__z^2)^(1/2))*(_i*k__y^2+_i*k__z^2-_j*k__x*k__y-_k*k__x*k__z)/((k__y^2+k__z^2)^(1/2)*((k__y^2+k__z^2)/k__y^2)^(1/2)*k__y)

(19)

Vector maths is really hard!

 

NULL

 

Download Perpendicular_Vectors_Rodrigues.mw

@ecterrab 

Magnificent!  Thank you Edgardo!

I am doing my best to understand the diferrences between the ChangeBasis and ChangeCoordinates commands.  Do you know where I can find help on the ChangeCoordinates command please?

Kind regards,

Steve

@Steve Roper Hi Edgardo,

I have found a way to acheive the result I was expecting, but I'm not really sure how it works.  This was inspired by your post here: https://www.mapleprimes.com/posts/212253-Vectors-In-Spherical-Coordinates-Using

Calculating dot products using spherical coordinates - Voodoo magic version````

with(Physics); with(Vectors)NULL

Setup(mathematicalnotation = true)

NULL

This is the first application of Voodoo (inspired by Edgardo's post here https://www.mapleprimes.com/posts/212253-Vectors-In-Spherical-Coordinates-Using)...

Setup(dimension = 3, coordinates = cartesian, g_ = `+`, spacetimeindices = lowercaselatin)

[coordinatesystems = {X}, dimension = 3, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}, spacetimeindices = lowercaselatin]

(1)

This is without voodoo...

k_ := _phi*c+_r*a+_theta*b

_phi*c+_r*a+_theta*b

(2)

``

k_.k_

a^2+b^2+c^2

(3)

I was hoping for the dot product to be calculated using the spherical coordinate system, but this is the cartesian dot product.NULL

 

Now for some more Voodoo...

tr := `~`[`=`]([X], ChangeCoordinates([X], spherical))

[x = r*sin(theta)*cos(phi), y = r*sin(theta)*sin(phi), z = r*cos(theta)]

(4)

cartesian_ := _i*x+_j*y+_k*z

_i*x+_j*y+_k*z

(5)

spherical_ := _phi*phi+_r*r+_theta*theta

_phi*phi+_r*r+_theta*theta

(6)

cartesian_ := subs(tr, cartesian_)

r*sin(theta)*cos(phi)*_i+r*sin(theta)*sin(phi)*_j+r*cos(theta)*_k

(7)

NULLNULL

cartesian_.cartesian_

r^2*sin(theta)^2*cos(phi)^2+r^2*sin(theta)^2*sin(phi)^2+r^2*cos(theta)^2

(8)

simplify(r^2*sin(theta)^2*cos(phi)^2+r^2*sin(theta)^2*sin(phi)^2+r^2*cos(theta)^2)

r^2

(9)

It works!

 

 

Download Calculating_dot_products_using_spherical_coordinates_with_Voodoo_majic.mw

 

 

@Carl Love Hi Carl,

Yes, &dot_sph is delivering the results that I was expecting from Maple's dot operator.

Thanks for your help.

1 2 Page 1 of 2