diff --git a/program2D/main.cc b/program2D/main.cc
index 125327575abe9e75a8909036c6e1058a9168fe78..d71342be32f14e3723491171ec7c499c8b3acb6b 100644
--- a/program2D/main.cc
+++ b/program2D/main.cc
@@ -17,23 +17,6 @@
using std::complex;
using namespace ::_COLSAMM_;
-/*
-std::complex<double> gauss(double x, double y) {
- //gauss electric fielöd
- double k = 2.0 * M_PI / lambda;
- double waist = radiusGauss * sqrt(1+pow(distanceFromWaist / rayleighrange,2));
- if (curvature == 0.0 || std::isnan(curvature))
- {
- return exp(-(x*x+y*y)/(waist*waist));
- }
- else
- {
-// return exp(-(x*x+y*y)/(waist*waist))*expi(-1.0 * k * (x*x+y*y) / (2.0 * curvature));
- }
-}
-*/
-
-
double fRealPart(std::complex< double > x) {
return x.real();
}
@@ -42,7 +25,11 @@ double fImagPart(std::complex< double > x) {
return x.imag();
}
-std::complex< double > fexpi(double x) {
+double absSquared(std::complex<double> value) {
+ return value.real()*value.real() + value.imag()*value.imag();
+}
+
+std::complex< double > expi(double x) {
std::complex< double > I(0.0,1.0);
return exp(I*x);
}
@@ -51,64 +38,200 @@ double fexp(double x) {
return exp(x);
}
+double mySqrt(std::complex<double> v) {
+ return sqrt(v.real());
+}
+
+void spectrumPlaneWave(double distance,
+ double wavenumber,
+ VariableFFT& varE,
+ VariableFFT& Kx,
+ VariableFFT& Ky) {
+ Blockgrid2D* blockGridRec = varE.Give_blockgrid();
+
+ X_coordinate2d X(*blockGridRec);
+ Y_coordinate2d Y(*blockGridRec);
+
+ double dx = varE.getHx();
+ double dy = varE.getHy();
+
+ double ukx = 2.0 * M_PI / varE.getSizeX();
+ double uky = 2.0 * M_PI / varE.getSizeY();
+
+ Kx = X / dx * ukx ;
+ Ky = Y / dy * uky ;
+
+ Function2d1<double, std::complex<double> > Sqrt(mySqrt);
+ Function2d1<std::complex<double>, double > Expi(expi);
+
+
+ varE.FFT();
+//Test NOW
+ varE = varE * Expi(-Sqrt(wavenumber*wavenumber - Kx*Kx - Ky*Ky) * distance);
+//varE = Expi(Sqrt(wavenumber*wavenumber - Kx*Kx - Ky*Ky) * distance);
+ varE.inversFFT();
+
+
+//varE = Sqrt(wavenumber*wavenumber - Kx*Kx - Ky*Ky);
+
+std::ofstream DATEIA;
+ DATEIA.open("varFarOut_test.vtk");
+ varE.Print_VTK(DATEIA);
+ DATEIA.close();
+
+
+}
+
+
+void CalcFresnel(double distance,
+ double lambda,
+ VariableFFT& varE) {
+ Blockgrid2D* blockGridRec = varE.Give_blockgrid();
+
+ int Nx = varE.getNx();
+ double sizeX = varE.getSizeX();
+
+ double wavenumber = 2.0 * M_PI / lambda;
+
+ X_coordinate2d X(*blockGridRec);
+ Y_coordinate2d Y(*blockGridRec);
+
+ Function2d1<std::complex<double>,double> Expi(expi);
+
+ varE = varE * Expi(-wavenumber * (X*X+Y*Y) / (2.0 * distance));
+ varE.FFT();
+
+ double strech = distance * lambda * Nx / (sizeX*sizeX);
+ varE.setDomainStrech(strech,strech);
+}
+
+
+
double SG_extern;
-double stretch_extern;
+double wSG_extern;
double rSG(double x, double y) {
double r = sqrt(x*x+y*y);
- r = r / stretch_extern;
+ r = r / wSG_extern;
return pow(r,SG_extern);
}
int main(int argc, char** argv) {
- std::ofstream DATEI;
+ cout.precision(10);
+ cout.setf(std::ios::fixed,std::ios::floatfield);
+ std::ifstream PARAMETER;
+ ofstream DATEI;
+
+ int levelFFT;
+ int typ; // 0: plane wave, 1: Fresnel
+ double zScreen;
+ double radiusGeo;
+
+
+ PARAMETER.open("para.dat",std::ios::in);
+ PARAMETER >> levelFFT >> typ >> zScreen >> radiusGeo;
+ PARAMETER.close();
+
+ // feste input parameter
double lambda = 1064e-6; //[mm]
+ double angle = 2.0; //[grad]
+ double zw = -2.0; //[mm]
+ double SG = 15.6;
+
+ // berechnete Werte
+ double angleRad = angle/180.0*M_PI;
+ double a = -zw;
+ double b = lambda / (M_PI * tan(angleRad) * tan(angleRad));
+ double R = (a*a+b*b)/a;
+ double wavenumber = 2.0 * M_PI / lambda;
+ double w = sqrt(2.0 * (a*a+b*b) / (b * wavenumber));
+ double wSG = sqrt(2.0) * w;
+
+ cout << " a: " << a << " b: " << b << endl
+ << " w: " << w << " wSG: " << wSG << endl
+ << " R: " << R
+ << endl;
+
+ SG_extern = SG;
+ wSG_extern = wSG;
+
Function2d1< double, std::complex< double > > realPart(fRealPart);
Function2d1< double, std::complex< double > > imagPart(fImagPart);
- Function2d1< std::complex< double >, double > Expi(fexpi);
+ Function2d1< std::complex< double >, double > Expi(expi);
Function2d1< double, double > Exp(fexp);
Function2d2< double, double > powSG(rSG);
+ Function2d1<double, std::complex<double> > AbsSquared(absSquared);
+
+ double h = radiusGeo / (pow(2,levelFFT));
+ Rechteck geo(-radiusGeo, -radiusGeo, radiusGeo-h, radiusGeo-h);
- double radiusGeo = 0.1;
-// double radiusGeo = 1.0;
- Rechteck geo(-radiusGeo, -radiusGeo, radiusGeo, radiusGeo);
-
- double w = 0.0705;
- stretch_extern = 1.41421;
-// double w = 0.5;
- double SG = 15.6;
-// double SG = 2.0;
- SG_extern = SG;
- double R = 2.03;
-
- double wavenumber = 2.0 * M_PI / lambda;
- int depth = 8;
-
- VariableFFT testU (depth, depth, geo);
- Blockgrid2D* bg = testU.Give_blockgrid();
+ VariableFFT varE (levelFFT, levelFFT, geo);
+ VariableFFT Kx(varE);
+ VariableFFT Ky(varE);
+
+ Blockgrid2D* bg = varE.Give_blockgrid();
+
+ Variable2D<double> powerField(*bg);
Variable2D<double> RealPart(*bg);
Variable2D<double> ImagPart(*bg);
X_coordinate2d X(*bg);
Y_coordinate2d Y(*bg);
- testU = Exp(-powSG(X/w,Y/w));
-// testU = Exp(-powSG(X,Y)/(w*w));
-
- testU = testU * Expi(-wavenumber*(X*X+Y*Y)/(2.0*R) );
+ //input
+ varE = Exp(-powSG(X,Y));
+ varE = varE * Expi(-wavenumber*(X*X+Y*Y)/(2.0*R) );
-// testU = Y;
+ //varE = X;
+
+ DATEI.open("varEin.vtk");
+ varE.Print_VTK(DATEI);
+ DATEI.close();
+
+ RealPart = realPart(varE);
+ ImagPart = imagPart(varE);
- DATEI.open("testU.vtk");
- testU.Print_VTK(DATEI);
+ DATEI.open("realUin.vtk");
+ RealPart.Print_VTK(DATEI);
DATEI.close();
- RealPart = realPart(testU);
- ImagPart = imagPart(testU);
+ DATEI.open("imagUin.vtk");
+ ImagPart.Print_VTK(DATEI);
+ DATEI.close();
+
+
+ if(typ==0) {
+ spectrumPlaneWave(zScreen,
+ wavenumber,
+ varE, Kx, Ky);
+ }
+ else {
+ CalcFresnel(zScreen,lambda,varE);
+ }
+
+ powerField = AbsSquared(varE);
+ ImagPart = varE.getHx() * varE.getHy();
+
+ double power = product(powerField,ImagPart);
+ powerField = powerField / power;
+
+ DATEI.open("powerField.vtk");
+ powerField.Print_VTK(DATEI);
+ DATEI.close();
+
+ /*
+ DATEI.open("varEout.vtk");
+ varE.Print_VTK(DATEI);
+ DATEI.close();
+ */
+
+
+ RealPart = realPart(varE);
+ ImagPart = imagPart(varE);
DATEI.open("realU.vtk");
RealPart.Print_VTK(DATEI);
@@ -117,6 +240,7 @@ int main(int argc, char** argv) {
DATEI.open("imagU.vtk");
ImagPart.Print_VTK(DATEI);
DATEI.close();
-
+
+
}