[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;
void STFTa(double f0,double t0,double w0,double x0,double k0,double w1,double k1,double s1x,double th,double xs1,double xd1);
void STFT(int Nt,double t[],double p1td[],double dt,double th,double t0,double w0);
void STFTn1(int Nt,double t[],double t0,double w0,double s1x,double xs1,double xd1,double k1,double w1,double th, double dt);
void imprime(int Nx,double x[],double C[],double p1x[]);
void grafica(int it);
int main()
{
double f0=1000,t0=1.0/f0,w0=2.0*pi*f0,x0=c*t0,k0=w0/c,alfa=1.3278;
double w1=2.0*w0,k1=w1/c,s1x=0.3*x0,th=0.3*t0;
int ix,Nx0=40,Nx=402,Ns1,Nd1;double dx=x0/double(Nx0),xs1,xd1,x[Nx+1],C[Nx+1],p1x[Nx+1];
int it,Nt=600; double dt=dx/(2.0*c),t[Nt+1],p1xd[Nt+1],p1td[Nt+1];
double aux1;
Ns1 = Nx/2; if ( Ns1 % 2 == 0)Ns1 = Ns1+1; xs1 = -double(Nx/2)*dx+double(Ns1)*dx;
xd1 = 3.0*x0; Nd1=int(xd1/dx + double(Nx)/2.0); if ( Nd1 % 2 == 0)Nd1 = Nd1+1;
for (ix=0 ; ix<=Nx ; ix++){x[ix]=-(double(Nx)/2.0)*dx+double(ix)*dx; C[ix]=0.0;}
for (it=0 ; it<=Nt ; it++){t[it]=-(double(Nt)/2.0)*dt+double(it)*dt;
p1td[it] = exp( -(0.5/(s1x*s1x))*(xd1-xs1-c*t[it])*(xd1-xs1-c*t[it]) )*cos(k1*(xd1-xs1)-w1*t[it]); // detector
}
STFTa(f0,t0,w0,x0,k0,w1,k1,s1x,th,xs1,xd1);
STFTn1(Nt,t,t0,w0,s1x,xs1,xd1,k1,w1,th,dt);
// -> 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]-xs1-c*t[it])*(x[ix+1]-xs1-c*t[it]) )*cos(k1*(x[ix+1]-xs1)-w1*t[it]);
if(ix+1==Nd1)
{
cout << "hola" << it+1 << "\n";
p1xd[it+1] = p1x[Nd1];
}
}
aux1 = exp( -(0.5/(s1x*s1x))*(xs1-xs1-c*t[it])*(xs1-xs1-c*t[it]) )*cos(k1*(xs1-xs1)-w1*t[it]);
C[Ns1] = C[Ns1] + aux1;
for (ix=0;ix<=Nx-4;ix+=2)
{
C[ix+2] = C[ix+2] – (0.5/alfa)*(C[ix+3]-C[ix+1]);
}
// <- fdtd
imprime(Nx,x,C,p1x);
grafica(it);
} //n
// <- time loop
ofstream myfile400;myfile400.open ("data400.dat");
for (int it=0;it<=Nt-2;it++)
{myfile400 << it << " " << t[it] << " " << p1td[it] << " " << p1xd[it] << "\n"; }
myfile400.close();
STFT(Nt,t,p1xd,dt,th,t0,w0);
}
////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////
void STFTa(double f0,double t0,double w0,double x0,double k0,double w1,double k1,double s1x,double th,double xs1,double xd1)
{
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);
double beta2_1,gamma_1,delta_1;
double aux1,aux2,aux3,aux4,aux5,FF1r,FF1i,FF1;
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 = ((xd1-xs1)*c)/(s1x*s1x)+tau/(th*th);
delta_1 = ((xd1-xs1)*(xd1-xs1))/(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*(xd1-xs1)+(gamma_1/beta2_1)*(w-w1);
aux5 = -k1*(xd1-xs1)+(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();
}
////////////////////////////////
////////////////////////////////
void STFTn1(int Nt,double t[],double t0,double w0,double s1x,double xs1,double xd1,double k1,double w1,double th,double dt)
{
int it;
double p1t0[Nt+1],p1td[Nt+1];
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);
double h,FF1r,FF1i,FF1;
double aux1,aux2;
ofstream myfile105;myfile105.open ("data105.dat");
for (it=0 ; it<=Nt ; it++)
{
p1t0[it] = exp( -(0.5/(s1x*s1x))*(xs1-xs1-c*t[it])*(xs1-xs1-c*t[it]) )*cos(k1*(xs1-xs1)-w1*t[it]); // origen
p1td[it] = exp( -(0.5/(s1x*s1x))*(xd1-xs1-c*t[it])*(xd1-xs1-c*t[it]) )*cos(k1*(xd1-xs1)-w1*t[it]); // detector
myfile105 << t[it] << " " << p1t0[it] << " " << p1td[it] << "\n";
}
myfile105.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))*((xd1-xs1)-c*t[it])*((xd1-xs1)-c*t[it]) );
aux2 = cos(k1*(xd1-xs1)-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();
}
////////////////////////////////
void imprime(int Nx,double x[],double C[],double p1x[])
{
ofstream myfile200;myfile200.open ("data200.dat");
for (int ix=0;ix<=Nx-2;ix+=2){myfile200 << x[ix+1] << " " << C[ix+1] << " " << p1x[ix+1] << "\n"; }
myfile200.close();
}
////////////////////////////////
void grafica(int it)
{
cout << it << "\n";
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");
}
//////////////////////////////////
////////////////////////////////
void STFT(int Nt,double t[],double p1td[],double dt,double th,double t0,double w0)
{
int it;
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);
double h,FF1r,FF1i,FF1;
double aux1,aux2;
ofstream myfile40;myfile40.open ("data160.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=1;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;
myfile40 << tau*1000.0 << " " << (w/(2.0*pi))*(1.0/1000.0) << " " << FF1 << "\n";
}
}
myfile40.close();
}
////////////////////////////////
[/code]