#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define KE 200
int main (int argc, const char * argv[])
{
float dx[KE], ex[KE], hy[KE],ix[KE];
float ga[KE], gb[KE];
int n,m,k,kc,kstart,nsteps;
float ddx,dt,T,epsz,epsilon,sigma;
float t0,spread,pi,pulse;
FILE *fp,*fopen();
float ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2;
float real_pt[5][KE],imag_pt[5][KE];
float freq[5],arg[5],ampn[5][KE],phasen[5][KE];
float real_in[5],imag_in[5],amp_in[5],phase_in[5];
float mag[KE];
kc=KE/2;
pi=3.14159;
epsz=8.8e-12;
ddx=0.01;
dt=ddx/(6e8);
printf("%6.4e %10.5e\n",ddx,dt);
for (k=1; k<KE; k++) {
ga[k]=1.;
gb[k]=0.;
dx[k]=0.;
ex[k]=0.;
hy[k]=0.;
ix[k]=0.;
mag[k]=0.;
// line 47
for (m=0; m<=2; m++) {
real_pt[m][k]=0.; // Real part and
imag_pt[m][k]=0.; // Imaginary part or the Fourier transform
ampn[m][k]=0.; // Amplitud y fase de las
phasen[m][k]=0.; // transformadas de Fourier.
}
}
for (m=0; m<=2; m++) {
real_in[m]=0.;
imag_in[m]=0.;
}
ex_low_m1=0.;
ex_low_m2=0.;
ex_high_m1=0.;
ex_high_m2=0.;
// Parámetros para la transformada de Fourier
freq[0]=100.0e6;
freq[1]=200.0e6;
freq[2]=500.0e6;
for (n=0; n<=2; n++) {
arg[n]=2*pi*freq[n]*dt;
printf("%2d %6.2f %7.5f \n",n,freq[n]*1e-6,arg[n]);
}
printf("Dielectric start at --->");
scanf("%d",&kstart);
printf("Epsilon --->");
scanf("%f",&epsilon);
printf("Conductivity --->");
scanf("%f",&sigma);
printf("%d %6.2f %6.2f \n",kstart,epsilon,sigma);
for (k=kstart; k<=KE; k++) {
ga[k]=1./(epsilon+sigma*dt/epsz);
gb[k]=sigma*dt/epsz;
}
for (k=1; k<=KE; k++) {
printf("%2d %6.2f %6.4f\n",k,ga[k],gb[k]);
}
t0=50.0;
spread=10.0;
T=0;
nsteps=1;
while (nsteps>0) {
printf("nsteps --->");
scanf("%d",&nsteps);
printf("%d \n",nsteps);
for (n=1; n<=nsteps; n++) { // Inicio del FOR principal
T=T+1;
// Calculo del campo Dx
for (k=0; k<KE; k++) {
dx[k]=dx[k]+0.5*(hy[k-1]-hy[k]);
}
// Inicializar con un pulso
pulse=exp(-0.5*(pow((t0-T)/spread, 2.0)));
dx[5]=dx[5]+pulse;
printf("%5.1f %6.2f %6.2f\n",T,pulse,dx[5]);
// Calcular el campo Ex a partir de Dx
for (k=0; k<KE-1; k++) {
ex[k]=ga[k]*(dx[k]-ix[k]);
ix[k]=ix[k]+gb[k]*ex[k];
}
// Calculo de la transformada de Fourier de Ex
for (k=0; k<KE; k++) {
for (m=0; m<=2; m++) {
real_pt[m][k]=real_pt[m][k]+cos(arg[m]*T)*ex[k];
imag_pt[m][k]=imag_pt[m][k]-sin(arg[m]*T)*ex[k];
}
}
// Transformada de Fourier del pulso de entrada
if (T<100) {
for (m=0; m<=2; m++) {
real_in[m]=real_in[m]+cos(arg[m]*T)*ex[10];
imag_in[m]=imag_in[m]-sin(arg[m]*T)*ex[10];
}
}
// Condiciones de frontera
ex[0] = ex_low_m2;
ex_low_m2 = ex_low_m1;
ex_low_m1 = ex[1];
ex[KE-1] = ex_high_m2;
ex_high_m2 = ex_high_m1;
ex_high_m1 = ex[KE-2];
// Calculo del campo Hy
for (k=0; k<KE-1; k++) {
hy[k]=hy[k]+0.5*(ex[k]-ex[k+1]);
}
}// Final del FOR principal
// Escribir el campo E a un archivo Ex
fp = fopen("Ex_f.dat","w");
for (k=1; k<=KE; k++) {
fprintf(fp, "%d %6.2f\n",k,ex[k]);
}
fclose(fp);
// Calcular la amplitud y la fase para cada frecuencia
// Amplitud y fase del pulso de entrada
for (m=0; m<=2; m++) {
amp_in[m]=sqrtf(pow(imag_in[m], 2.))+pow(real_in[m], 2.);
phase_in[m]=atan2(imag_in[m], real_in[m]);
printf("%d Input pulse : %8.4f %8.4f %8.4f %7.2f\n",m,real_in[m],imag_in[m],amp_in[m],(180.0/pi)*phase_in[m]);
for (k=1; k<KE; k++) {
ampn[m][k]=(1./amp_in[m])*sqrt(pow(real_pt[m][k], 2.)+pow(imag_pt[m][k], 2.));
phasen[m][k]=atan2(imag_pt[m][k], real_pt[m][k])-phase_in[m];
}
// Escribir la amplitud del campo al archivo Amp
fp=fopen("Amp0.dat", "w");
for (k=0; k<KE; k++) {
fprintf(fp, " %d %8.5f",k,ampn[0][k]);
}
fclose(fp);
fp=fopen("Amp1.dat", "w");
for (k=0; k<KE; k++) {
fprintf(fp," %d %8.5f \n",k,ampn[1][k]);
}
fclose(fp);
fp=fopen("Amp2.dat", "w");
for (k=0; k<KE; k++) {
fprintf(fp," %d %8.5f \n",k,ampn[2][k]);
}
fclose(fp);
printf(" %5.1f\n",T);
}
} // Final del LOOP While
}
评论0