first commit

Maxwell3D.m

@@ -0,0 +1,82 @@
+%physical constants
+clear all;
+c    = 2.998e8;
+eta0 = 120*pi;
+mu0  = pi*4e-7;
+eps0 = 1e-9/(36*pi);
+%environment parameters
+width  = 1;
+height = 1;
+depth  = 1;
+f0     =1e9;
+tw     = 1e-8/pi;%?
+t0     = 4*tw;
+%discretization parameters
+dx = 0.01;
+dy = dx;
+dz = dx;
+nx = width/dx;
+ny = height/dy;
+nz = depth/dz;
+dt   = 0.95/(c*sqrt(dx^-2+dy^-2+dz^-2));
+%calculation parameters
+n_iter = 2000;
+%initalization
+Hx = zeros(nx,ny,nz);
+Hy = zeros(nx,ny,nz);
+Hz = zeros(nx,ny,nz);
+Ex = zeros(nx,ny,nz);
+Ey = zeros(nx,ny,nz);
+Ez = zeros(nx,ny,nz);
+%iteration
+i = 0
+for n=1:1:n_iter
+    %derivatives
+
+    %H
+    Hxy = diff(Hx,1,2);
+    Hxz = diff(Hx,1,3);
+    Hzx = diff(Hz,1,1);
+    Hzy = diff(Hz,1,2);
+    Hyx = diff(Hy,1,1);
+    Hyz = diff(Hy,1,3);
+    %Maxwell Equations
+    % boyutlar e? uzunlukta olmas? için türev al?nmayan boyutu bir eksik
+    % indeksliyoruz. E vektörlerinin türev boyutunu k?saltmadan al?p di?rev
+    % vektörün türev boyutunda bir ileriden ba?lat?yoruz.
+
+    %E
+    Ex(:,2:end-1,2:end-1) = Ex(:,2:end-1,2:end-1) + (dt/(eps0*dy))*Hzy(:,1:end-1,2:end-1) - (dt/(eps0*dz))*Hyz(:,2:ny-1,1:end-1); %OK
+    Ey(2:end-1,:,2:end-1) = Ey(2:end-1,:,2:end-1) + (dt/(eps0*dz))*Hxz(2:end-1,:,1:end-1) - (dt/(eps0*dx))*Hzx(1:end-1,:,2:end-1); %OK
+    Ez(2:end-1,2:end-1,:) = Ez(2:end-1,2:end-1,:) + (dt/(eps0*dx))*Hyx(1:end-1,2:end-1,:) - (dt/(eps0*dy))*Hxy(2:end-1,1:end-1,:); %OK
+     %Gaussian Source
+      f(n)= sin(2*pi*f0*n*dt)*exp(-(n*dt-t0)^2/(tw^2))/dy;
+    Ez(round(nx/2),round(ny/2),round(nz/2)) = Ez(round(nx/2),round(ny/2),round(nz/2)) + f(n);
+    Ezn(n)=Ez(round(nx/2),round(ny/2)-4,round(nz/2));
+    Exy = diff(Ex,1,2);
+    Exz = diff(Ex,1,3);
+    Ezx = diff(Ez,1,1);
+    Ezy = diff(Ez,1,2);
+    Eyx = diff(Ey,1,1);
+    Eyz = diff(Ey,1,3);
+
+    %H
+    Hx(:,1:end-1,1:end-1) = Hx(:,1:end-1,1:end-1) - (dt/(mu0*dy))*Ezy(:,:,1:end-1) + (dt/(mu0*dz))*Eyz(:,1:end-1,:); %OK
+    Hy(1:end-1,:,1:end-1) = Hy(1:end-1,:,1:end-1) - (dt/(mu0*dz))*Exz(1:end-1,:,:) + (dt/(mu0*dx))*Ezx(:,:,1:end-1); %OK
+    Hz(1:end-1,1:end-1,:) = Hz(1:end-1,1:end-1,:) - (dt/(mu0*dx))*Eyx(:,1:end-1,:) + (dt/(mu0*dy))*Exy(1:end-1,:,:); %OK
+    %display
+    if (mod(i,5)==0)
+        aa(:,:)= Ez(:,round(ny/2),:);
+        a = log(aa.*(aa>0));
+        a(1,1)= -60;
+        a(1,2)=0;
+        a(end,1)=-60;
+        a(end,end)= -60;
+        a(1,end)=-60;
+       pcolor(a);colorbar
+    end
+    %shading interp;
+    drawnow
+    i = i+1
+end
+

README.md

@@ -0,0 +1 @@
+# TR-FTDT

bitirme1.m

@@ -0,0 +1,52 @@
+%TMz Polarization
+%physical constants
+c    = 2.998e8;
+eta0 = 120*pi;
+mu0  = pi*4e-7;
+eps0 = 1e-9/(36*pi);
+%environment parameters
+width  = 1;
+height = 1;
+tw     = 1e-9/pi;%?
+t0     = 4*tw;
+%discretization parameters
+dx = 0.005;
+dy = 0.005;
+nx = width/dx;
+ny = height/dy;
+dt   = 0.95/(c*sqrt(dx^-2+dy^-2));
+%calculation parameters
+n_iter = 1000;
+c1 = dt/(mu0*dy);
+c2 = dt/(mu0*dx);
+c3 = dt/(eps0*dx);
+c4 = dt/(eps0*dy);
+%initalization
+Hx = zeros(nx,ny);
+Hy = zeros(nx,ny);
+Ez = zeros(nx,ny);
+%iteration
+for n=1:1:n_iter
+    %Maxwell Equations (TMz)
+    Ezx = diff(Ez,1,1);
+    Ezy = diff(Ez,1,2);
+    Hx(2:nx-1,2:ny) = Hx(2:nx-1,2:ny) - c1*Ezy(2:nx-1,:);
+    Hy(2:nx,2:ny-1) = Hy(2:nx,2:ny-1) + c2*Ezx(:,2:ny-1);
+    Hxy = diff(Hx,1,2);
+    Hyx = diff(Hy,1,1);
+    Ez(2:nx-1,2:ny-1) = Ez(2:nx-1,2:ny-1) + c3*Hyx(2:nx-1,2:ny-1) - c4*Hxy(2:nx-1,2:ny-1);
+    %Gaussian Source
+    f(n)= exp(-(n*dt-t0)^2/(tw^2))/dy;
+    Ez(round(nx/2),round(ny/2)) = Ez(round(nx/2),round(ny/2)) + f(n);
+    %Neuman Condition
+    Ez(:,2)  = -Ez(:,1);
+    Ez(2,:)  = -Ez(1,:);
+    Ez(:,ny-1) = -Ez(:,ny);
+    Ez(nx-1,:) = -Ez(nx,:);
+    %display
+    %n = n + 1;
+    pcolor(Ez);
+    shading interp;
+    drawnow
+end
+

maxwell3d.py

@@ -0,0 +1,103 @@
+import torch
+from math import sin, exp, sqrt, pi
+
+def get_device():
+    try:
+        device = torch.device('cuda')
+        assert device
+    except:
+        device = torch.device('cpu')
+    return torch.device('cpu') #device
+
+def zeros(m,n,k):
+    return torch.zeros((m,n,k), device=get_device(), dtype=torch.double)
+
+def diff(a,dim):
+    if dim == 0:
+        return a[1:,:,:] - a[:-1,:,:]
+    elif dim == 1:
+        return a[:,1:,:] - a[:,:-1,:]
+    elif dim == 2:
+        return a[:,:,1:] - a[:,:,:-1]
+
+#physical constants
+c    = 2.998e8
+eta0 = 120*pi
+mu0  = pi*4e-7
+eps0 = 1e-9/(36*pi)
+#environment parameters
+width  = 1
+height = 1
+depth  = 1
+f0     = 1e9
+tw     = 1e-8/pi
+t0     = 4*tw
+#discretization parameters
+dx = 0.01
+dy = 0.01
+dz = 0.01
+nx = width  / dx
+ny = height / dy
+nz = depth  / dz
+dt   = 0.95 / (c*sqrt(dx**-2+dy**-2+dz**-2))
+#calculation parameters
+n_iter = 100
+#initalization
+Hx = zeros(nx,ny,nz)
+Hy = zeros(nx,ny,nz)
+Hz = zeros(nx,ny,nz)
+Ex = zeros(nx,ny,nz)
+Ey = zeros(nx,ny,nz)
+Ez = zeros(nx,ny,nz)
+#iteration
+i = 0
+
+for n in range(n_iter):
+
+    #derivatives
+
+    Hxy = diff(Hx,1);
+    Hxz = diff(Hx,2);
+    Hzx = diff(Hz,0);
+    Hzy = diff(Hz,1);
+    Hyx = diff(Hy,0);
+    Hyz = diff(Hy,2);
+
+    #Maxwell Equations  
+
+    Ex[:,1:-1,1:-1] += (dt/(eps0*dy))*Hzy[:,:-1,1:-1] - (dt/(eps0*dz))*Hyz[:,1:-1,:-1]
+    Ey[1:-1,:,1:-1] += (dt/(eps0*dz))*Hxz[1:-1,:,:-1] - (dt/(eps0*dx))*Hzx[:-1,:,1:-1]
+    Ez[1:-1,1:-1,:] += (dt/(eps0*dx))*Hyx[:-1,1:-1,:] - (dt/(eps0*dy))*Hxy[1:-1,:-1,:]
+
+    #Gaussian Source
+
+    Ez[int(nx/2),int(ny/2),int(nz/2)] += sin(2*pi*f0*n*dt)*exp(-(n*dt-t0)**2/(tw**2))/dy
+
+    #derivatives
+
+    Exy = diff(Ex,1);
+    Exz = diff(Ex,2);############################!!!!!!!!!!!
+    Ezx = diff(Ez,0);
+    Ezy = diff(Ez,1);
+    Eyx = diff(Ey,0);
+    Eyz = diff(Ey,2);
+
+    #Maxwell Equations 
+
+    Hx[:,:-1,:-1] += -(dt/(mu0*dy))*Ezy[:,:,:-1] + (dt/(mu0*dz))*Eyz[:,:-1,:]
+    Hy[:-1,:,:-1] += -(dt/(mu0*dz))*Exz[:-1,:,:] + (dt/(mu0*dx))*Ezx[:,:,:-1]
+    Hz[:-1,:-1,:] += -(dt/(mu0*dx))*Eyx[:,:-1,:] + (dt/(mu0*dy))*Exy[:-1,:,:]
+
+    ####################### I AM HERE #########################
+
+    """
+    #display
+    if not i%5:
+       pcolor(Ez[:,int(ny/2),:])
+    drawnow
+    """
+
+    i += 1
+    print(i)
+
+