Skip to main content

FDTD Method

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

Download original notebook

Define 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[*)(Graphics[Line[Offload[electricField]], PlotRange -> {{0, 400}, {-1, 1}}, TransitionDuration -> 5])(*,*)(*"1:eJxdTVsKAjEQqy/UD3+8gXcR9WNBKV6g7k41UFrpdO/hjTVFZHHnI5NkmGR3S9ZPjDE6J5xS6PysqhXhmN3zgVaHe4MoX7UknL0PyXW6IZcgbcloDxBGTH8Ptg+ia5JLSMW6eJfh1kDLvwK3wYsFI//NQe0dJW9JrtlFRUGK+z67urGg/QE7azJt"*)(*]VB*)

Manipulate the time

EventHandler[InputRange[1,800,10,10], Function[val, 
  electricField = {Range[400], generator[val]} // Transpose;
]]
(*VB[*)(EventObject[<|"Id" -> "c89daae1-37e8-42c0-8b1d-df9245dee993", "Initial" -> 10, "View" -> "8f7a33e3-c23f-4773-9462-7d69117199b2"|>])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKW6SZJxobpxrrJhsZp+mamJsb61qamBnpmqeYWRoamhtaWiYZAQB7JhTJ"*)(*]VB*)
Last updated on