fft.tar.gz_fft

/* Fd2d_eigen.c 2D program to find TM eigenfrequencies using fft */ # include <math.h> # include <stdlib.h> # include <stdio.h> # include <string.h> #define IE 102 #define JE 62 int main () { int l, n, i, j, f, ic, jc, nsteps = 0, Nhan, isource, jsource, T = 0, size; double eps[IE][JE]; double han[30000]; double ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE]; double real_pt[IE][JE], imag_pt[IE][JE], amp[IE][JE], phase[IE][JE]; double real_pt_fft[IE][JE], imag_pt_fft[IE][JE], amp_fft[IE][JE], phase_fft[IE][JE], amp_fft_rl[IE][JE]; double ddx, dt, epsz, pi, epsilon, sigma, eaf; double t0, spread, pulse; double xdist, ydist, dist; double stepfc, fc_step, sfc, fc_s, efc, fc_e, fc, arg; FILE *fp, *fp_fft, *fp2_fft, *fp_fft_rl; ic = IE/2; jc = JE/2; ddx = .001; dt = ddx/6e8; epsz = 8.8e-12; pi=3.14159; for ( j=0; j < JE; j++ ) { for ( i=0; i < IE; i++ ) { dz[i][j] = 0.; ez[i][j] = 0.; hx[i][j] = 0.; hy[i][j] = 0.; ga[i][j]= 1.0 ; eps[i][j] = 1.; real_pt[i][j] = 0.; imag_pt[i][j] = 0.; real_pt_fft[i][j] = 0.; imag_pt_fft[i][j] = 0.; amp_fft[i][j] = 0.; phase_fft[i][j] = 0.; amp_fft_rl[i][j] = 0.; } } /* Specify the dielectric in the waveguide*/ for ( j=0; j < 40; j++ ) { for ( i=0; i < 30; i++ ) { eps[i][j] = 3.0; ga[i][j] = 1./eps[i][j]; } } for ( i=60; i < IE; i++ ) { for ( j= 40; j < JE; j++ ) { eps[i][j] = 2.0; ga[i][j] = 1./eps[i][j]; } } /* Initialize the test function- Gaussian Waveform */ spread = 7.0; for ( j=1; j < JE; j++ ) { for ( i=1; i < IE; i++ ) { isource = ic+10-i; jsource = jc-10-j; xdist = isource; ydist = jsource; dist = sqrt(pow(xdist,2.) + pow(ydist,2.)); dz[i][j] = exp(-.5*(pow((dist/spread),2.))); } } isource = ic + 9; jsource = jc - 11; printf( "nsteps -- > ");//20000 scanf("%d", &nsteps); // Use Nhan = 20000 to find the eigenfrequencies and Nhan = 5000 to find the eigenfunctions. Nhan = nsteps; for ( n=0; n < Nhan; n++ ) han[n] = .5*(1 - cos(2*pi*n/(Nhan-1))); printf( "start fc (GHz) -- > ");//0 scanf("%lf", &sfc); fc_s = sfc*1e9; printf( "end fc (GHz) -- > ");//6 scanf("%lf", &efc); fc_e = efc*1e9; printf( "step fc (GHz) -- > ");//0.003 scanf("%lf", &stepfc); fc_step = stepfc*1e9; fp = fopen( "Ez_time","w"); for ( n=0; n <=nsteps ; n++) { /*----Start of the Main FDTD loop ----*/ /* Calculate the Dz field */ for ( j=1; j < JE; j++ ) for ( i=1; i < IE; i++ ) dz[i][j] = dz[i][j]+ .5*( hy[i][j] - hy[i-1][j] - hx[i][j]+ hx[i][j-1]) ; /* Calculate the Ez field */ /* The field is truncated to simulate metal boundaries */ for ( j=2; j < JE-1; j++ ) for ( i=2; i < IE-1; i++ ) ez[i][j] = ga[i][j]*dz[i][j] ; /* Calculate the Hx field */ for ( j=0; j < JE-1; j++ ) for ( i=0; i < IE-1; i++ ) hx[i][j] = hx[i][j] + .5*( ez[i][j] - ez[i][j+1] ) ; /* Calculate the Hy field */ for ( j=0; j < JE-1; j++ ) for ( i=0; i <= IE-1; i++ ) hy[i][j] = hy[i][j] + .5*( ez[i+1][j] - ez[i][j] ) ; // Time Domain data at source position to find eigenfrequencies using fourier transform fprintf(fp,"%lf\n", ez[isource][jsource]); T = T + 1; }/*----End of the main FDTD loop ---- */ fclose(fp); fp = fopen( "Time","w"); fprintf(fp,"T = %10.5f ps \n",1e9*dt*(float)T); fclose(fp); // Reading Time domain data to memory int lines_allocated = 128; int max_line_len = 100; /* Allocate lines of text */ char **ezr = (char **)malloc(sizeof(char*)*lines_allocated); fp = fopen("Ez_time", "r"); int m; for (m=0;1;m++) { /* Have we gone over our line allocation? */ if (m >= lines_allocated) { int new_size; /* Double our allocation and re-allocate */ new_size = lines_allocated*2; ezr = (char **)realloc(ezr,sizeof(char*)*new_size); lines_allocated = new_size; } /* Allocate space for the next line */ ezr[m] = malloc(max_line_len); if (fgets(ezr[m],max_line_len-1,fp)==NULL) break; /* Get rid of CR or LF at end of line */ int p; for (p=strlen(ezr[m])-1; p>=0 && (ezr[m][p]=='\n' || ezr[m][p]=='\r'); p--); ezr[m][p+1]='\0'; } fclose(fp); //Fast fourier transform calculations fc = fc_s; size = ceil((fc_e - fc_s)/fc_step); T = 0; fp_fft = fopen( "Amp_fft","w"); fp2_fft = fopen( "Phase_fft","w"); fp_fft_rl = fopen( "Amp_fft_real","w"); for( f = 0; f <= size; f++) { for ( n=0; n <=nsteps ; n++) { arg = 2*pi*dt*fc; real_pt_fft[isource][jsource] = real_pt_fft[isource][jsource] + han[n]* cos(arg*T) * strtod(ezr[n], NULL); imag_pt_fft[isource][jsource] = imag_pt_fft[isource][jsource] - han[n]* sin(arg*T) * strtod(ezr[n], NULL); T = T + 1; } amp_fft[isource][jsource] = sqrt( pow(real_pt_fft[isource][jsource],2.) + pow(imag_pt_fft[isource][jsource],2.) ); phase_fft[isource][jsource] = atan2(imag_pt_fft[isource][jsource],real_pt_fft[isource][jsource]); amp_fft_rl[isource][jsource] = amp_fft[isource][jsource]*cos(phase_fft[isource][jsource]); fprintf( fp_fft,"%lf\t%lf\n",fc, amp_fft[isource][jsource]); fprintf( fp2_fft,"%lf\t%lf\n",fc, phase_fft[isource][jsource]); fprintf( fp_fft_rl,"%lf\t%lf\n",fc, amp_fft_rl[isource][jsource]); fc = fc + fc_step; } fclose(fp_fft); fclose(fp2_fft); fclose(fp_fft_rl); for (;m>=0;m--) free(ezr[m]); free(ezr); return 0; }