Need help for FDTD code for wave propagation !!!!!!!!!!! I am unable to detect errors in my code.

This code is for gaussian wave propagating in step index fiber(with core and cladding) This is 2D FDTD simulation code. the wave is travelling along z and the Gaussian beam is input source(varying along y). TE mode of propagation is considered. Please tell any error in the code.(Why the light is not confined in the core region, )

clc;
warning off all;   

%*********************************************************
%      Fumdamental constants
%**********************************************************
  c=2.99792458 * 10^8 * 10^6;  %speed of light
  mu0=4*(22/7)*10^-7 * 10^-6;  %permeability
  epsilon0=1/(c.^2 * mu0);  %permiitivity
  FSI=sqrt(mu0/epsilon0);    %free space impedance
  lambda=1.55;      %wavlength

%**********************************************************
%      Structure dimensions
%***********************************************************
  Y=2;      %transverse length is structure
  Z=1;      %length of waveguide
  wow=.3;      %width of waveguide
  dz=lambda/200;
  dy=lambda/200;
  dx=lambda/200;
  dt0=1/(c*(1/dy^2 + 1/dz^2)^0.5);  %time resolution
  dt=dt0*.9;        %dt must less than dt0 for numerical stability condtion to be satisfied
  %**********************Derived Dimensions********************
  y=[-(Y/2-mod(Y/2,dy)):dy:(Y/2)];  %used for plotting
  z=[0:dz:Z];        %used for plotting

%*********************************************************
%      Material properties
%*********************************************************
  n=ones(length(z)+2,length(y)+2);  %refractive index profile
    %create the refractive index profile here
    %
    for aa=1:length(y);
      if(abs(y(aa)-dy) <= (wow/2) )
        n(:,aa)=3.5;   %core refractive index
      else
        n(:,aa)=1.5;   %cladding refractive index
      end
    end
    
    
    %
    %
  epsilon=(n.^2)*epsilon0;    %derived permitivity profile
  mu=ones(length(z)+2,length(y)+2)*mu0;  %premeability profile
  
  S=0*ones(length(z)+2,length(y)+2);  %electrical conductivity profile
  R=0*ones(length(z)+2,length(y)+2);  %magnetic resistance profile

%************************************************************
%      Field initialisation to 0
%**************************************************************

  Ex=zeros(length(z)+2,length(y)+2);  Exs=size(Ex);    
  Ey=zeros(length(z)+2,length(y)+2);  Eys=size(Ey);
  Hz=zeros(length(z)+2,length(y)+2);  Hzs=size(Hz);

  
ctr=10;  %just used as counter
for zzz=1:50000     %to run the loop for some time
    %**************************************************************
    %        Source
    %**************************************************************

          %gaussian wave input
          Hz(1,:)=[0,1*exp(-(y/(wow/exp(1))).^2   ),0];

  %************************************************************
  %      Field calculations
  %********************************************************************

    for aa=2:Exs(1,1)-1;
      for bb=2:Exs(1,2)-1;
        Ex(aa,bb)= (dt/epsilon(aa,bb)) * ( Hz(aa,bb)-Hz(aa,bb-1) )/dy + ( 1-S(aa,bb)*dt/epsilon(aa,bb) ) * Ex(aa,bb);
        Ey(aa,bb)=-(dt/epsilon(aa,bb)) * ( Hz(aa,bb)-Hz(aa-1,bb) )/dx + ( 1-S(aa,bb)*dt/epsilon(aa,bb) ) * Ey(aa,bb);
      end
    end
    

    for aa=2:Hzs(1,1)-1;
      for bb=2:Hzs(1,2)-1;
        Hz(aa,bb)= (dt/mu(aa,bb)) * ( .5*(Ex(aa,bb+1) -Ex(aa,bb))/dy -.5*(Ey(aa+1,bb)-Ey(aa,bb))/dx ) + (1-R(aa,bb)*dt/mu(aa,bb)) * Hz(aa,bb);
      end
    end

  %*******************************************************
  %      PLOTTING THE GRAPH
  %*******************************************************
    view(90,90);
    figure(1);
    mesh(y,z,Hz(2:Hzs(1,1)-1,2:Hzs(1,2)-1));
  
    
    ctr=ctr-1;
    if(ctr<=0 )
      ctr=input('Continue for how much time ??');
      
      if(ctr<=0)
        break;
      end
    end
     
end

The main problem in the above code is the the light is not confined in the core region. the V parameter is satisfied for the single mode. Still then, as the wave travels, the intensity of wave goes on decreasing.

This code is needed in my ongoing project for waveguide simulation. I m stuck with the problem for a long time. I might have to leave this project if i am not able to detect the error .

Please HELP, its urgent.

Answer 1

[Saurabh Suman]

code can also be downloaded from http://www.filefactory.com/file/ceef3be/n/waveguide.m

[Jan]

Most users will not download your code. It is your part of the question versus answer design to post all relevant information in a format, which is as easy to read as possible.

Answer 2

[Jan]

warning off all is a really bad idea. Warnings help to locate bugs and you are looking for bugs. As long as you ignore warnings, it does not seem, like this problem is "urgent" for your.

I do not know, what "the light is not confined in the core region" might mean. It is your turn to explain the problem such, that we do not have to start some investigations about the same topic.

I see your program and you claim, that there is a problem with the light. But there is no variable called "light" or anything equivalent. Therefore I have no chance to understand what the problem is and in consequence I cannot suggest a solution.

Please take the time to explain yor problem again by editing the original message. Makeing your code more readable is a good idea also.

[Saurabh Suman]

Thanks for the reply. And i do accept your point.

Now i m checking with warning on.

The above code is complete and coded based on Finite Difference Time Domain method for solving Maxwell's Equations, obtained from Yee's algorithm. Ex, Ey and Hz are 3 light components. And according to code, confinement of Hz component is being looked after.

[adi]

hi, im looking for a very similar code to what you wrote and i can see what you are saying about the light beeing non confined in the waveguide core. did you have any solution for this problem?