/*
 *                     WAVES IN TWO DIMENSIONS
 */
#include"lib351.h"

void main(){

    // program variables initialization
    		
    ofstream dataOutStr;
    outFileMaker("m353p1c.out", &dataOutStr);
    
    const double M_PI = float(3.14159);
    const double xScale = 0.78;
    const double yScale = 0.78;
    const double thres[7] = {0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32};
    const char pchar[7] =  {' ', 'I', '-', 'X', 'O', '#', 'M'};
	const int imColumns = 87;
    const int imRows = 55;
    int flag, i, j, n, index, run;
    double disc, l1, l2, s1, s2, v, r1, r2, d, ratio;
    double delta, x0;
    double cons, a, b, c, a1, a2, x, y, ai, aj;
	double rint[imColumns];
	char image[imColumns*13];

	// cycle through the projects
  	bool goout = false;
	run = 0;
	while (!goout){

    	run++;

		// data input

    	cout<<"Please, supply values :"<<endl;
    	cout<<"\nseparation between sources, delta ="<<endl;
    	cin>>delta;
    	cout<<"\nratio of two source strenghts, ratio ="<<endl;
    	cin>>ratio;
    	cout<<"\nwavelengths"<<endl;
    	cout<<"\nl1 ="<<endl;
    	cin>>l1;
    	cout<<"\nl2"<<endl;
    	cin>>l2;
    	cout<<"\nvelocity of the source one, v ="<<endl;
    	cin>>v;

        // print the header
        
		if (run > 1) {
			dataOutStr<<"\f";
		}
		dataOutStr<<"WAVES IN TWO DIMENSIONS"<<endl;
        dataOutStr<<"delta = "<<delta<<";   ratio = "<<ratio;
        dataOutStr<<"; l1="<<l1<<"; l2="<<l2<<"; v="<<v<<endl;
    	dataOutStr<<"xScale ="<<xScale<<"; yScale ="<<yScale<<endl<<endl;

    	// for each row: fill in the columns and print the row

    	s1 = 1;
    	s2 = ratio;

    	for(i = 0; i < imRows; i++){
      		for(j = 0; j < imColumns; j++){
				ai = i;
				aj = j;
				x = aj /(5*xScale) - 51./5;
				y = ai /(3*yScale) - 31./3;

				r2 = sqrt(pow(x, 2) + pow(y-delta/2 ,2));
				if(v == 0){
	  				r1 = sqrt(pow(x, 2)+pow(y+delta/2 ,2));
				}
				else{
	    			d = 3 * v;
	  				if(d > 10) {
	    				d = 10.0;
	    			}
        			a = 1 - 1./(v*v);
        			b = 2*d/(v*v) - 2*x;
        			c = x*x + y*y - pow((d/v), 2);
        			if (v == 1) {
        				x0 = -c/b;
        			}
        			else {
        				disc = b*b - 4*a*c;
        				r1 = l1;
        				if(disc >= 0){
        					x0 = (-b + sqrt(disc))/(2*a);
        					if(x0 <= d){
	    						r1 = sqrt(pow((x - x0), 2) + y*y);
	    					}
	  					}
	  				}
				}

				r1 += 0.000001;
				r2 += 0.000001;
				l2 += 0.000001;

				a1 = s1 * sin(2*M_PI*r1/l1)/r1;
				a2 = s2 * sin(2*M_PI*r2/l2)/r2;

				cons = 2./(s1*s1 + s2*s2);
				rint[j] = cons * pow((a1+a2), 2);

				if(v == 0){
	  				if (r1 <= 0.7){
	    				rint[j] = 0;
	  				}
	  				if (s2 != 0 && r2 <= 0.7){
	    				rint[j] = 0;
	  				}
            	}
      		}

        	// digitize array 'rint' as array 'image'

      		for(j = 0; j < imColumns*13; j++){
				image[j] = ' ';
      		}
      		index = 0;
      		for(n = 0; n < 7; n++){
				flag = 0;
     			for(j = 0; j < imColumns; j++){
	  				index++;
	  				if(rint[j] >= thres[n]){
	    				image[index] = pchar[n];
	  				}
				}
				if (n<=6) {
					for (j = 0; j < imColumns; j++) {
						index++;
						image[index] = '\b';
					}
				}
      		}
  			for (j = 0; j < imColumns*13 ; j++){
          		dataOutStr<<image[j];
        	}
  			dataOutStr<<endl;
    	}
    	goout = query();
  	}
    fileClose("m353p1c.out", &dataOutStr);
};