Solving Maxwell's equation for 1D case with two interfaces. Based on the work of John B. Schneider "Understanding the Finite-Difference Time-Domain Method" 2023 we will try to solve 1D Maxwell's equation for electro-magnetic wave propagating through a non-dispersive plate
Download original notebookDefine constants
ζ = 377.0; (*BB[*)(*Ohm*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) size = 400; ez = ConstantArray[0., {size}]; hy = ConstantArray[0., {size}]; cμ = ConstantArray[ζ, {size}]; cϵ = ConstantArray[1., {size}]; (*BB[*)(*Define a second medium with losses*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) cϵ[[200;;220]] = (*FB[*)(((1 - 0.001))(*,*)/(*,*)((1 + 0.001)))(*]FB*); cμ[[200;;220]] = (*FB[*)((ζ)(*,*)/(*,*)(9.0))(*]FB*) / (1 + 0.001); generator = With[{ie = size}, Compile[{{steps, _Integer}}, Module[ {ez = ez, hy = hy}, Do[ (*BB[*)(*Calculate magnetic field *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) Do[ hy[[j]] += (*FB[*)(((ez[[j + 1]] - ez[[j]]))(*,*)/(*,*)(ζ))(*]FB*); , {j, 1, ie-1}]; (*BB[*)(*Boundary conditions*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) hy[[ie]] = hy[[ie-1]]; ez[[ie]] = ez[[ie-1]]; (*BB[*)(*Pulse generator*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) hy[[49]] -= Exp[- (*FB[*)(((*SpB[*)Power[(n - 30.)(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)(100.))(*]FB*)] / ζ; ez[[1]] = ez[[2]]; (*BB[*)(*Calculate electric field *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) Do[ ez[[j]] = cϵ[[j]] ez[[j]] + cμ[[j]] (hy[[j]] - hy[[j - 1]]) , {j, 2, ie - 1}]; (*BB[*)(*Pulse generator*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) ez[[50]] += Exp[- (*FB[*)(((*SpB[*)Power[(n + 0.5 - (-0.5) - 30.)(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)(100.0))(*]FB*)]; , {n, steps}]; ez], CompilationOptions -> {"InlineExternalDefinitions" -> True}, "CompilationTarget" -> "C", "RuntimeOptions" -> "Speed"]];
Plot as a function of x-coordinate
electricField = {Range[400], generator[3]} // Transpose; Graphics[ Line[electricField // Offload], PlotRange->{{0,400}, {-1,1}}, TransitionDuration->5 ]
(*VB[*)(FrontEndRef["c163d4ea-9002-4a25-a9f0-9e37f0a858d1"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxuaGaeYpCbqWhoYGOmaJBqZ6iZaphnoWqYam6cZJFqYWqQYAgCBDhVu"*)(*]VB*)
Try to drag this slider 😉
Manipulate the time
EventHandler[InputRange[1,800,20,1], Function[val, electricField = {Range[400], generator[val]} // Transpose; ]]
(*VB[*)(EventObject[<|"Id" -> "97462fa3-8bc0-472a-919b-54cc650eacf1", "Initial" -> 1, "View" -> "f294fa8d-e0cd-416c-9be1-baffe32749b6"|>])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKpxlZmqQlWqTophokp+iaGJol61ompRrqJiWmpaUaG5mbWCaZAQCVOxZr"*)(*]VB*)
Congratulations! 🥂 Now you know how to solve 1D Maxwell's equation with an arbitary boundary conditions!