Simulacion

 

 

[code language=”cpp”]
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <complex>
#include <time.h>
using namespace std;

double pi = 2.0*asin(1.0);
double c = 343.0;

int main()
{
double f0=1000,t0=1.0/f0,w0=2.0*pi*f0,x0=c*t0,k0=w0/c,th=0.3*t0,alfa=1.3278;
int ix,Nx=402,Nx0=40,Ns1,Nd1;

double dx=x0/double(Nx0),x[Nx+1],C[Nx+1];

int it,Nt=600;
double dt=dx/(2.0*c),t[Nt+1],h;

double p1tO[Nt+1],p1td[Nt+1],p1x[Nx+1],s1x=0.3*x0,w1=2.0*w0,k1=w1/c,x1,xd,p1xd[Nt+1];

double beta2_1,gamma_1,delta_1,FF1r,FF1i,FF1;
double aux1,aux2,aux3,aux4,aux5;
int itau,Ntau=50;double tau,taui=-1.0*t0,tauf=+4.0*t0,dtau=(tauf-taui)/float(Ntau);
int iw, Nw =50;double w, wi =+0.0*w0,wf =+5.0*w0,dw =(wf-wi) /float(Nw);

Ns1 = Nx/2;
if ( Ns1 % 2 == 0)Ns1 = Ns1+1;
x1 = -double(Nx/2)*dx+double(Ns1)*dx;

xd = 3.0*x0;
tau = 0.0*t0;

Nd1=int(xd/dx + double(Nx)/2.0);
if ( Nd1 % 2 == 0)Nd1 = Nd1+1;

cout << "La extension de la malla es " << double(Nx)*dx << "\n";
cout << "Limite inf " << -double(Nx/2)*dx << " y sup " << +double(Nx/2)*dx << "\n";
cout << "La malla tiene " << Nx << " particiones " << "\n";
cout << "Los limites son " << -double(Nx)/2 << " y " << double(Nx)/2 << "\n";
cout << "La fuente esta en " << x1 << " en la particion " << Ns1 << "\n";
cout << "\n";
cout << "El detector esta en " << Nd1 << "\n";
cout << "\n";
cout << "frecuencia f0 = "<< f0 << " " << "\n";
cout << "tiempo referencia t0=1/f0 = "<< t0 << " " << "\n";
cout << "distancia de referencia x0=c*t0 = " << x0 << "\n";
cout << "dx =" << dx << "\n";
cout << "dt =" << dt << "\n";
cout << "dt del fdtd " << x0/(20.0*2.0*c) << "\n";
cout << "\n";
cout << "Datos del pulso " << "\n";
cout << "s1x " << s1x << "\n";
cout << "x1 " << x1 << "\n";
cout << "f1 " << w1/(2.0*pi) << "\n";
cout << "k1 " << k1 << "\n";

ofstream myfile100;myfile100.open ("data100.dat");
for (ix=0 ; ix<=Nx ; ix++)
{
x[ix] = -(double(Nx)/2.0)*dx + double(ix)*dx;
C[ix] = 0.0;
p1x[ix] = exp( -(0.5/(s1x*s1x))*(x[ix]-x1)*(x[ix]-x1) )*cos(k1*(x[ix]-x1));
myfile100 << x[ix] << " " << p1x[ix] << "\n";
}
myfile100.close();

ofstream myfile105;myfile105.open ("data105.dat");
for (it=0 ; it<=Nt ; it++)
{
t[it] = -(double(Nt)/2.0)*dt + double(it)*dt;
p1tO[it] = exp( -(0.5/(s1x*s1x))*(x1-x1-c*t[it])*(x1-x1-c*t[it]) )*cos(k1*(x1-x1)-w1*t[it]); // origen
p1td[it] = exp( -(0.5/(s1x*s1x))*(xd-x1-c*t[it])*(xd-x1-c*t[it]) )*cos(k1*(xd-x1)-w1*t[it]); // detector
h = exp( -(0.5/(th*th)) *(t[it]-tau)*(t[it]-tau) );
p1xd[it] = 0.0;
myfile105 << t[it] << " " << p1tO[it] << " " << p1td[it] << " " << h << "\n";
}
myfile105.close();

// -> time loop
for (it=0 ; it<= Nt ; it+=2)
{
// -> fdtd
for (ix=0;ix<=Nx-2;ix+=2)
{
C[ix+1] = C[ix+1] – 0.5*alfa*(C[ix+2]-C[ix]);
p1x[ix+1] = exp( -(0.5/(s1x*s1x))*(x[ix+1]-x1-c*t[it])*(x[ix+1]-x1-c*t[it]) )*cos(k1*(x[ix+1]-x1)-w1*t[it]);
if(ix+1==Nd1)
{
cout << "hola" << it+1 << "\n";
p1xd[it+1] = p1x[Nd1];
}
}
C[Ns1] = C[Ns1] + p1tO[it];
for (ix=0;ix<=Nx-4;ix+=2)
{
C[ix+2] = C[ix+2] – (0.5/alfa)*(C[ix+3]-C[ix+1]);
}
// <- fdtd

// -> output data
ofstream myfile200;myfile200.open ("data200.dat");
for (ix=0;ix<=Nx-2;ix+=2){myfile200 << x[ix+1] << " " << C[ix+1] << " " << p1x[ix+1] << "\n"; }
// for (ix=0;ix<=Nx-2;ix+=2){myfile200 << x[ix+1] << " " << p1x[ix+1] << "\n"; }
myfile200.close();

// ofstream myfile201;myfile201.open ("data201.dat");
// for (ix=0;ix<=Nx;ix+=2) {myfile201 << x[ix] << " " << C[ix] << "\n";}
// myfile201.close();
// <- output data

// -> output grafico
stringstream ss;ss << it;string fileName00="";fileName00="plot00.plt";
ofstream outputFile00;
outputFile00.open(fileName00.c_str());
outputFile00 << "set term png" << "\n";
outputFile00 << "set nokey" << "\n";
if(it<10) {outputFile00 << "set output \’a00"+ss.str()+".png\’" << "\n";}
if(it>=10 && it < 100) {outputFile00 << "set output \’a0"+ss.str()+".png\’" << "\n";}
if(it>=100 && it < 1000) {outputFile00 << "set output \’a"+ss.str()+".png\’" << "\n";}
outputFile00 << "set yra[-2:2]" << "\n";
outputFile00 << "set grid" << "\n";
outputFile00 << "plot ‘data200.dat’ u 1:2 w l,” u 1:3 w lp pt 7 ps 0.4 lc 2" << "\n";
// outputFile00 << "plot ‘data200.dat’ u 1:2 w lp" << "\n";
outputFile00.close();
system("gnuplot plot00.plt \n");
// <- output grafico
} //n

// <- time loop

system("convert a*.png a.gif \n");

ofstream myfile111;myfile111.open ("data111.dat");
for(itau=0;itau<=Ntau;itau++)
{
tau = taui+float(itau)*dtau;
beta2_1 = (c*c)/(s1x*s1x)+1.0/(th*th);
gamma_1 = ((xd-x1)*c)/(s1x*s1x)+tau/(th*th);
delta_1 = ((xd-x1)*(xd-x1))/(s1x*s1x)+(tau*tau)/(th*th);
for(iw=0;iw<=Nw;iw++)
{
w = wi+float(iw)*dw;
aux1 = exp( -0.5*delta_1+0.5*(gamma_1*gamma_1)/beta2_1 );
aux2 = exp( -0.5*((w-w1)*(w-w1))/beta2_1 );
aux3 = exp( -0.5*((w+w1)*(w+w1))/beta2_1 );
aux4 = +k1*(xd-x1)+(gamma_1/beta2_1)*(w-w1);
aux5 = -k1*(xd-x1)+(gamma_1/beta2_1)*(w+w1);
FF1r = (0.5/(sqrt(beta2_1)*th))*(aux1*(aux2*cos(aux4)+aux3*cos(aux5)));
FF1i = (0.5/(sqrt(beta2_1)*th))*(aux1*(aux2*sin(aux4)+aux3*sin(aux5)));
FF1 = FF1r*FF1r+FF1i*FF1i;
myfile111 << tau*1000.0 << " " << (w/(2.0*pi))*(1.0/1000.0) << " " << FF1 << "\n";
}
}
myfile111.close();

ofstream myfile30;myfile30.open ("data30.dat");
for(itau=0;itau<=Ntau;itau++)
{
tau = taui+float(itau)*dtau;
for(iw=0;iw<=Nw;iw++)
{
w = wi+float(iw)*dw;

FF1r = 0.0;
FF1i = 0.0;
for(it=0;it<=Nt;it++)
{
aux1 = exp( -(1.0/(2.0*s1x*s1x))*((xd-x1)-c*t[it])*((xd-x1)-c*t[it]) );
aux2 = cos(k1*(xd-x1)-w1*t[it]);
h = (1.0/(th*sqrt(2.0*pi)))*exp( -(0.5/(th*th))*(t[it]-tau)*(t[it]-tau) );
FF1r = FF1r + h*aux1*aux2*cos(w*t[it]);
FF1i = FF1i + h*aux1*aux2*sin(w*t[it]);
}
FF1r = dt*FF1r;
FF1i = dt*FF1i;
FF1 = FF1r*FF1r+FF1i*FF1i;
myfile30 << tau*1000.0 << " " << (w/(2.0*pi))*(1.0/1000.0) << " " << FF1 << "\n";
}
}
myfile30.close();

ofstream myfile40;myfile40.open ("data40.dat");
for(itau=0;itau<=Ntau;itau++)
{
tau = taui+float(itau)*dtau;
for(iw=0;iw<=Nw;iw++)
{
w = wi+float(iw)*dw;

FF1r = 0.0;
FF1i = 0.0;
for(it=0;it<=Nt;it++)
{
h = (1.0/(th*sqrt(2.0*pi)))*exp( -(0.5/(th*th))*(t[it]-tau)*(t[it]-tau) );
FF1r = FF1r + h*p1td[it]*cos(w*t[it]);
FF1i = FF1i + h*p1td[it]*sin(w*t[it]);
}
FF1r = dt*FF1r;
FF1i = dt*FF1i;
FF1 = FF1r*FF1r+FF1i*FF1i;
myfile40 << tau*1000.0 << " " << (w/(2.0*pi))*(1.0/1000.0) << " " << FF1 << "\n";
}
}
myfile40.close();

ofstream myfile50;myfile50.open ("data50.dat");
for(itau=0;itau<=Ntau;itau++)
{
tau = taui+float(itau)*dtau;
for(iw=0;iw<=Nw;iw++)
{
w = wi+float(iw)*dw;

FF1r = 0.0;
FF1i = 0.0;
for(it=0 ; it<= Nt ; it+=2)
{
h = (1.0/(th*sqrt(2.0*pi)))*exp( -(0.5/(th*th))*(t[it]-tau)*(t[it]-tau) );
FF1r = FF1r + h*p1td[it]*cos(w*t[it]);
FF1i = FF1i + h*p1td[it]*sin(w*t[it]);
}
FF1r = 2.0*dt*FF1r;
FF1i = 2.0*dt*FF1i;
FF1 = FF1r*FF1r+FF1i*FF1i;
myfile50 << tau*1000.0 << " " << (w/(2.0*pi))*(1.0/1000.0) << " " << FF1 << "\n";
}
}
myfile50.close();
//

ofstream myfile205;myfile205.open ("data205.dat");
for (it=0 ; it<=Nt ; it++)
{
myfile205 << t[it] << " " << p1tO[it] << " " << p1td[it] << " " << p1xd[it] << "\n";
}
myfile205.close();

//

}

[/code]

 

Aqui esta una version con mejor estilo en C++

d