Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Need help with my 2D FDTD code written in C language

Status
Not open for further replies.

Miriam Griffin Hera

Newbie level 2
Newbie level 2
Joined
Aug 21, 2014
Messages
2
Helped
0
Reputation
0
Reaction score
0
Trophy points
1
Visit site
Activity points
9
Hi, can anyone help me with my 2D FDTD code that I wrote using C language? Any help will be greatly appreciated.

Below is the code:

Code:
# include <math.h>
# include <stdlib.h>
# include <stdio.h>

# define IE 60
# define JE 60

int main ()
{
  float dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
  int l, n, i, j, ic, jc, nsteps;
  float dx, dy, dt, T, epsz, pi, sigma, mu;
  float t0, spread, pulse;
  FILE *fp,  *fopen();

  ic  =  IE/2;
  jc  =  JE/2;
  dx  =  .01;        /* cell size */
  dy  =  .01;        /* cell size */
  dt  =  dx/6e8;     /* Time steps */
  epsz  =  8.8e-12;
  pi  = 3.14159;
  mu  = 1.0e-7*4.0*pi;
  sigma  = 1.0e-8;

  for( j = 0; j < JE; j++) {
    printf("%2d", j);
    for(i = 0; i < IE; i++) {
      ez[i][j]  =  0.;
      hx[i][j]  =  0.;
      hy[i][j]  =  0.;
    }
    printf("\n");
  }

  t0 = 20.0;
  spread = 6.0;
  T = 0;
  nsteps = 1;

  while(nsteps > 0){
    printf("nsteps -> ");
    scanf("%d", &nsteps);
    printf("%d \n", nsteps);

    for(n = 1; n <= nsteps; n++) {
      T  = T;

        /*  - Start of Main FDTD loop  - */

      /*  -  Put a Gaussian pulse in the middle  -  */
    
      pulse = sin(2*pi*1500*1e6*dt*T);
      ez[ic][jc] = pulse;

      /*  -  Calculate Ez field  -  */

     for(j = 1; j < JE; j++){
        for(i = 1; i < IE; i++) {
          ez[i][j] = (1.0 - sigma*dt/(2.0*epsz))/(1.0+sigma*dt)/(2.0*epsz)*ez[i][j]+(dt/epsz)
            /(1.0+sigma*dt/(2.0*epsz*(1/dx))*(hy[i][j] - hy[i - 1][j]) - (dt/epsz)/((1.0+sigma*dt)/(2.0*epsz))*(1/dy)
            *(hx[i][j] - hx[i][j - 1]));
        }
      }

      /*  -  Calculate Hx field  -  */
      for(j = 0; j < JE - 1; j++){
        for(i = 0; i < IE - 1; i++){
          hx[i][j] = hx[i][j] - (dt/mu*1/dy)*(ez[i][j+1] - ez[i][j]);
        }
      }
      
      /* Calculate Hy field */
      for(j = 0; j < JE - 1; j++){
        for(i = 0; i <= IE - 1; i++) {
          hy[i][j] = hy[i][j]+(dt/mu*1/dx)*(ez[i+1][j] - ez[i][j]);
        }
      }
    }

    /*  -  End of main FDTD loop  -  */

    for(j = 1; j < jc; j++) {
      printf("%2d", j);
      for(i = 1; i < ic; i++) {
        printf("%5.2f", ez[2*i][2*j]);
      }
      printf("\n");
    }
    printf("T = %5.0f \n", T);

    /* Write E field out to file "Ez" */

    fp = fopen("Ez", "w");
    for(j = 0; j < JE; j++) {
      for(i = 0; i < IE; i++) {
        fprintf(fp,"%d\t%d\t%6.3f\n",i, j, ez[i][j]);
      }
      fprintf(fp, "\n");
    }
    fclose(fp);
  }
}

Thank you!
 

Status
Not open for further replies.

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top