// interferenz.cc
// 
// Berechnung von Interferenzbildern

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <complex>
#include <vector>
#include <string>
#include <sstream>
#include <algorithm>
#include <climits>
#include "mpi.h"

using namespace std;

const double two_pi=6.283185307179586;

// Struktur zur Beschreibung eines Bildschirms
typedef struct {
  double x1, y1, x2, y2;
} screen;

// Struktur zur Beschreibung der Aufloesung des Bildschirms
typedef struct {
  int x, y;
} resolution;

void read_input_file(const string &filename, double &dist, screen &S,
		     resolution &res, double &lambda, 
		     vector<double> &x, vector<double> &y, 
		     vector<double> &amplitude) {
  // sehr einfacher Parser
  ifstream in(filename.c_str());
  do {
    string token;
    in >> token;
    if (in) {
      if (token=="dist")
	in >> dist;
      else if (token=="screen") 
	in >> S.x1 >> S.y1 >> S.x2 >> S.y2;
      else if (token=="picture")
	in >> res.x >> res.y;
      else if (token=="lambda")
	in >> lambda;
      else if (token=="source") {
	double t_x, t_y, t_amplitude;
	in >> t_x >> t_y >> t_amplitude;
	x.push_back(t_x);
	y.push_back(t_y);
	amplitude.push_back(t_amplitude);
      };
      in.ignore(INT_MAX, '\n');
      // wenn #include <limits> besser
      // in.ignore(numeric_limits<int>::max(), '\n');
    }
  } while (in);
}

void write_output_file(const string &filename, resolution res,
		       const vector<double> &intensity,
		       const double max_intensity) {
  // schreibe Daten in ein pgm-File (portable graymap file format), hellster
  // Punkt auf Schirm wird weiss gefaerbt andere grau je nach Helligkeit
  ofstream out(filename.c_str());
  out << "P2" << endl  // magic number
      << res.x << " " << res.y << endl  // Breite und Hoehe 
      << 256 << endl;  // Anzahl der Grauwerte
  for (int i=0; i<res.x*res.y; ++i) 
    out << floor(intensity[i]/max_intensity*255.0+0.5) << endl;
}

void calc_intensity(const double dist, const screen S, const resolution res,
		    const double lambda, const vector<double> x, 
		    const vector<double> y, const vector<double> amplitude, 
		    vector<double> &intensity) {
  // brechnet fuer jeden Punkt (x, y) des Bildschirms die Intensitaet
  intensity.resize(res.x*res.y);
  double dx=S.x2-S.x1;
  double dy=S.y2-S.y1;
  double dist2=dist*dist;
  for (int j=0; j<res.y; ++j) {
    double Y=S.y1+(j+0.5)*dy/res.y;
    for (int i=0; i<res.x; ++i) {
      double X=S.x1+(i+0.5)*dx/res.x;
      complex<double> comp_int=(0.0, 0.0);
      for (int k=0; k<x.size(); ++k) { // Summation ueber alle Quellen
	double d=sqrt((X-x[k])*(X-x[k])+(Y-y[k])*(Y-y[k])+dist2);
	comp_int+=amplitude[k]*
	  exp(complex<double>(0.0, two_pi*d/lambda))/d;
      }
      // Intensitaet ist das Betragsquadrat der complexen Amplitude
      intensity[j*res.x+i]=norm(comp_int);
    }
  }
}

int main(void) {

  MPI::Init();  // initialisiere MPI

  int rank=MPI::COMM_WORLD.Get_rank();  // eigenen Rang und 
  int size=MPI::COMM_WORLD.Get_size();  // Zahl der Prozesse ermitteln

  // Parameter mit Defaultwerten belegen
  double dist=1.0;  // Abstand von Quellen zum Schirm  
  screen S={-0.05, -0.05,   // linke untere Ecke des Schirms
	    +0.05, +0.05};  // rechte obere Ecke des Schirms
  resolution res={250, 250};   // Anzahl der Punke in x- und y-Richtung
  double lambda=0.4e-6;  // Wellenlaenge der Strahlung
  vector<double> x, y, amplitude;
  int num_sources;

  if (rank==0) {  // Rang 0 liest Konfigurationsdatei ein
    if (argc==2) 
      read_input_file(argv[1], dist, S, res, lambda, x, y, amplitude);
    num_sources=x.size();
    if (num_sources==0)  // ohne Quellen ist nichts zuberechnen
      MPI::COMM_WORLD.Abort(EXIT_FAILURE);
  }
  
  // Prozess 0 macht Parameter allen anderen Prozessen bekannt
  MPI::COMM_WORLD.Bcast(&dist, 1, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&S.x1, 1, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&S.y1, 1, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&S.x2, 1, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&S.y2, 1, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&res.x, 1, MPI::INT, 0);
  MPI::COMM_WORLD.Bcast(&res.y, 1, MPI::INT, 0);
  MPI::COMM_WORLD.Bcast(&lambda, 1, MPI::DOUBLE, 0);    
  MPI::COMM_WORLD.Bcast(&num_sources, 1, MPI::INT, 0);
  // alle Prozesse ausser Prozess 0 muessen erst Speicher reservieren
  if (rank!=0) {
    x.resize(num_sources);
    y.resize(num_sources);
    amplitude.resize(num_sources);
  }
  MPI::COMM_WORLD.Bcast(&x[0], num_sources, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&y[0], num_sources, MPI::DOUBLE, 0);
  MPI::COMM_WORLD.Bcast(&amplitude[0], num_sources, MPI::DOUBLE, 0);

  // y-Aufloesung aufrunden, so dass duch Anzahl der Prozesse teilbar
  if (res.y%size>0)
    res.y=res.y-res.y%size+size;

  // von jedem Prozess zu berechnende Aufloesung
  resolution res_local={res.x, res.y/size};

  // von jedem Prozess zu berechnende Flaeche
  screen S_local={S.x1, rank*(S.y2-S.y1)/size+S.y1,
		  S.x2, (rank+1)*(S.y2-S.y1)/size+S.y1};
  
  // hier startet die Berechnung
  vector<double> intensity_local;
  calc_intensity(dist, S_local, res_local, lambda, x, y, amplitude, 
		 intensity_local);

  // maximale Intensitaet ermitteln, erst lokal dan global
  double max_intensity, max_intensity_local=
    *max_element(intensity_local.begin(), intensity_local.end());
  MPI::COMM_WORLD.Reduce(&max_intensity_local, &max_intensity, 1, MPI::DOUBLE, 
			 MPI::MAX, 0);

  // Prozess 0 sammelt alle Ergebnisse ein ...
  vector<double> intensity;
  intensity.resize(res.x*res.y);
  MPI::COMM_WORLD.Gather(&intensity_local[0], res_local.x*res_local.y, 
			 MPI::DOUBLE,
			 &intensity[0], res_local.x*res_local.y,
			 MPI::DOUBLE, 0);
  // ... und schreibt sie in eine Datei
  if (rank==0)
    write_output_file(string(argv[1])+".pgm", res, intensity, max_intensity);
  
  MPI::Finalize();  // beende MPI

  return EXIT_SUCCESS;
}