#include "declarations.h"
//#define _INTERACTIVE_
//#define _TEXOUTPUT_
#define _GNUPLOTOUTPUT_
double laplace2d(int x, int y, double** potential);
bool checkError(double oldPotential, double potential, double accuracy);
int main (int argc, char * const argv[]) {
//Set gnuplot path, comment out at home
Gnuplot::set_GNUPlotPath("/opt/local/bin");
//Open Files for editing
intermediateOut.open(intermediateFileName.c_str());
approxSolutionOut.open(approxSolutionFileName.c_str());
exactSolutionOut.open(exactSolutionFileName.c_str());
gnuplotOut.open(gnuplotFileName.c_str());
//Set output style
intermediateOut.setf(ios::showpoint|ios::fixed);
approxSolutionOut.setf(ios::showpoint|ios::fixed);
exactSolutionOut.setf(ios::showpoint|ios::fixed);
gnuplotOut.setf(ios::showpoint|ios::fixed);
//Set output file precision
intermediateOut.precision(2);
approxSolutionOut.precision(2);
exactSolutionOut.precision(2);
gnuplotOut.precision(2);
bool goodToGo = false;
//gnuplot objects for plots
Gnuplot surfacePlot("lines");
Gnuplot contourPlot("lines");
Gnuplot exactContourPlot("lines");
//seed random number
srand((unsigned)time(0));
#ifdef _INTERACTIVE_
//Collect paramaters from the user
cout << "Please enter a resolution: ";
cin >> resolution;
cout << endl << "Please enter an accuracy: ";
cin >> accuracy;
#endif
int arraySize = resolution + 2;
cout << endl << "Initilazing Potential Field." << endl;
//Initializing data arrays
//Final relaxation solution array
double** potential=new double*[arraySize];
for(i=0; i < arraySize; i++){
potential[i]=new double[arraySize];
for(j=0; j<arraySize; j++){
potential[i][j] = lowest+int(range*rand()/(RAND_MAX + 1.0));
}
}
//Intermediate relaxation temp array
double** oldPotential=new double*[arraySize];
for(i=0; i<arraySize; i++){
oldPotential[i]=new double[arraySize];
for(j=0; j<arraySize; j++){
oldPotential[i][j] = 0;
}
}
//Exact solution array
double** exactPotential=new double*[arraySize];
for(i=0; i<arraySize; i++){
exactPotential[i]=new double[arraySize];
for(j=0; j<arraySize; j++){
oldPotential[i][j] = 0;
}
}
//Set boundry conditions
for(i=0; i < arraySize; i++){
potential[16][i] = 100.00;
oldPotential[16][i] = 100.00;
}
for(j = 1; j < arraySize; j++){
potential[j][0] = 0.00;
oldPotential[j][0] = 0.00;
}
for(i = 0; i < arraySize; i++){
potential[0][i] = 0.00;
oldPotential[0][i] = 0.00;
}
for(j = 0; j < arraySize; j++){
potential[j][16] = 0.00;
oldPotential[j][16] = 0.00;
}
//data arrays ready to go
//calculate potential using laplaces equation
while(!goodToGo){
for(i = 1; i<= resolution; i++){
for(j = 1; j <= resolution; j++){
oldPotential[i][j] = potential[i][j];
potential[i][j] = laplace2d(i,j,potential);
}
}
//Determine if another iteration is needed
for(i = 1; i < resolution; i++){
for(j = 1; j < resolution; j++){
goodToGo = checkError(oldPotential[i][j], potential[i][j], accuracy);
if(!goodToGo){
break;
}
}
if(!goodToGo){
break;
}
}
}
//Calculate exact solution (first 225 terms)
double x = 0;
double y = 0;
for(i = 0; i < arraySize; i++){
for(j = 0; j < arraySize; j++){
for(int k = 0; k < 225; k++){
if( 0 == ((k+1) % 2)){
x = j/(double)(arraySize - 1);
y = i/(double)(arraySize - 1);
exactPotential[i][j] += (1/sinh(k*PI)) * sin((k*x*PI)) * sinh((k*y*PI))/k;
}
}
exactPotential[i][j] *= ((4*Vnot)/PI);
}
}
//output data to gnuplot data file
for(i = 0; i < arraySize; i++){
for(j=0; j < arraySize; j++){
gnuplotOut << potential[i][j] << " ";
}
gnuplotOut << endl;
}
#ifdef _GNUPLOTOUTPUT_
//Output for plotting
//Output exact solution
for(i = 0; i < arraySize; i++){
for(j = 0; j < arraySize; j++){
//Simple space delimited output.
exactSolutionOut << exactPotential[i][j] << " ";
}
exactSolutionOut << endl;
}
//Output approx solution
for(i = 0; i < arraySize; i++){
for(j = 0; j < arraySize; j++){
//Simple space delimited output.
approxSolutionOut << potential[i][j] << " ";
}
approxSolutionOut << endl;
}
//Output intermediate solution
for(i=resolution-1; i>=0; i--){
for(j=resolution-1; j>=0; j--){
//Simple space delimited output.
intermediateOut << oldPotential[i][j] << " ";
}
intermediateOut << endl;
}
#endif
//TeX Output -- Important! The indexs are switched around here
#ifdef _TEXOUTPUT_
//Output exact solution
for(i = 16; i >= 0; i--){
for(j = 0; j < arraySize; j++){
//TeX tabular output, plotting on this data will fail.
if(0 == j) exactSolutionOut << exactPotential[i][j];
else exactSolutionOut << exactPotential[i][j] << " & ";
}
//TeX tabular output, plotting on this data will fail.
exactSolutionOut << " \\\\ \\hline" << endl;
}
//Output approx solution
for(i = 16; i >= 0; i--){
for(j = 0; j < arraySize; j++){
//TeX tabular output, plotting on this data will fail.
if(0 == j) approxSolutionOut << potential[i][j];
else approxSolutionOut << potential[i][j] << " & ";
}
//TeX tabular output, plotting on this data will fail.
approxSolutionOut << " \\\\ \\hline" << endl;
}
//Output intermediate solution
for(i = 16; i >= 0; i--){
for(j=resolution-1; j>=0; j--){
//TeX tabular output, plotting on this data will fail.
if(0 == j) intermediateOut << oldPotential[i][j];
else intermediateOut << oldPotential[i][j] << " & ";
}
//TeX tabular output, plotting on this data will fail.
intermediateOut << " \\\\ \\hline" << endl;
}
#endif
#ifdef _GNUPLOTOUTPUT_
//This plot should be commented out if the latex output is turned on
//Plot exact solution contour plot, using gnuplot exceptions
try{
exactContourPlot.cmd("set ticslevel 0");
//exactContourPlot.cmd("set view 70, 120");
exactContourPlot.cmd("set zlabel \"V(x,y)\" 5,1");
exactContourPlot.set_xlabel("x");
exactContourPlot.set_ylabel("y");
//exactContourPlot.set_zlabel("\"V(x,y)\" 0,1");
exactContourPlot.set_xrange(0.0, 16.00);
exactContourPlot.set_yrange(0.0, 16.00);
exactContourPlot.cmd(("splot \"" + exactSolutionFileName + "\" matrix title \"Exact Solution\" with lines"));
}
catch (GnuplotException ge){
cout << ge.what() << endl;
}
#endif
//Plot surface plot, using gnuplot exceptions
try{
surfacePlot.cmd("set ticslevel 0");
surfacePlot.cmd("set zlabel \"V(x,y)\" 5,1");
surfacePlot.set_xlabel("x");
surfacePlot.set_ylabel("y");
surfacePlot.set_xrange(0.0, 16.00);
surfacePlot.set_yrange(0.0, 16.00);
surfacePlot.cmd(("splot \"" + gnuplotFileName + "\" matrix title \"Potential\" with lines"));
}
catch (GnuplotException ge){
cout << ge.what() << endl;
}
//Plot equipotential lines, using gnuplot exceptions
try{
contourPlot.set_contour("base");
contourPlot.cmd("set title \"Equipotential Lines\"");
contourPlot.cmd("set view map");
contourPlot.cmd("unset surface");
contourPlot.cmd("set cntrparam bspline");
contourPlot.cmd("set cntrparam levels disc 40,70,90");
contourPlot.cmd("set zlabel \"V(x,y)\" 5,1");
contourPlot.set_xlabel("x");
contourPlot.set_ylabel("y");
contourPlot.set_xrange(0.0, 16.00);
contourPlot.set_yrange(0.0, 16.00);
contourPlot.cmd(("splot \"" + gnuplotFileName + "\" matrix with lines"));
}
catch (GnuplotException ge){
cout << ge.what() << endl;
}
cout << endl << "That's all folks" << endl;
return 0;
}
//function to approximate the solution to the laplace equation around a point x,y
double laplace2d(int x, int y, double** potential){
return(.25*(potential[x+1][y]+potential[x-1][y]+potential[x][y+1]+potential[x][y-1]));
}
//Function to check the error between each iteration
bool checkError(double oldPotential, double potential, double accuracy){
if (abs((oldPotential - potential))/oldPotential <= accuracy){
return true;
}
else{
return false;
}
}