#include "iostream"
#include "fstream"
#include "iomanip"
#include "cmath"
#include "stdio.h"
#include "string.h"

#define G0 6.67191e-3
#define pi (4.0*atan(1.0))
#define MAX 1e+30
#define MIN -1e+30

#define PARA "-p"
#define POSI "-l"
#define RANGE "-r"
#define SPACE "-i"
#define OUTFILE "-o"

#define BOLDRED "\033[1m\033[31m" /* Bold Red */
#define RESET "\033[0m" /* reset */

using namespace std;

class vec_cy
{
public:
	vec_cy();
	~vec_cy();
	int info_taker(double*,double*,double*,double*,char*);
	double func(double);
	double integration(double,double,double);
	int calculate();
	int datout();
private:
	double t_dep,b_dep,rad,den;
	double x_center,y_center;
	double xmin,xmax,ymin,ymax;
	double x_interval,y_interval,inte_interval;
	char* outfile_name;

	int M,N;
	double* data;

	double cal_radius;
};

vec_cy::vec_cy()
{
	data = NULL;
	inte_interval = 1e-5;
	outfile_name = NULL;
}

vec_cy::~vec_cy()
{
	if(data!=NULL) delete []data;
}

int vec_cy::info_taker(double* para,double* posi,double* range,double* interval,char* filename)
{
	t_dep = *para; b_dep = *(para+1); rad = *(para+2); den = *(para+3);
	x_center = *posi; y_center = *(posi+1);
	xmin = *range; xmax = *(range+1); ymin = *(range+2); ymax = *(range+3);
	x_interval = *interval; y_interval = *(interval+1);
	if (*(interval+2)>0)	
		inte_interval = *(interval+2);
	outfile_name = filename;

	M = int (ymax-ymin)/y_interval+1;
	N = int (xmax-xmin)/x_interval+1;

	ymax = ymin + (M-1)*y_interval;
	xmax = xmin + (N-1)*x_interval;

	calculate();
	datout();
	return 0;
}

double vec_cy::func(double l)
{
	double temp,temp1;
	temp=1.0/sqrt(l*l+t_dep*t_dep)-1.0/sqrt(l*l+b_dep*b_dep);
  	temp1=2*l*acos((l*l+cal_radius*cal_radius-rad*rad)/(2*l*cal_radius));
  	if (cal_radius==rad) return (2*l*acos(l/(2*rad)))*temp;
  	else return temp*temp1;
}

double vec_cy::integration(double a,double b,double eps)
{
	int n,k;
	double fa,fb,h,t1,p,s,x,t;
	fa=func (a);  fb=func (b);
	n=1; h=b-a;
	t1=h*(fa+fb)/2.0;
	p=eps+1.0;
	while (p>=eps)
	{ 
		s=0.0;
		for (k=0;k<=n-1;k++)
		{ 
			x=a+(k+0.5)*h;
			s=s+func (x);
		}
		t=(t1+h*s)/2.0;
		p=fabs(t1-t);
		t1=t; n=n+n; h=h/2.0;
	}
	return t;
}

int vec_cy::calculate()
{
	double temp_res;
	data = new double [M*N];
	for (int i = 0; i < M*N; i++)
	{
		cal_radius = (int) sqrt(pow(xmin+x_interval*(i%N)-x_center,2)+pow(ymin+y_interval*(i/N)-y_center,2));
		if (cal_radius >= rad)
		{
			data[i] = G0*den*integration(cal_radius-rad,cal_radius+rad,inte_interval);
		}
		else if (cal_radius>0&&cal_radius<rad)
		{
			temp_res = 2*pi*(sqrt(pow(rad-cal_radius,2)+t_dep*t_dep)-sqrt(pow(rad-cal_radius,2)+b_dep*b_dep)+b_dep-t_dep);
			data[i] = G0*den*(temp_res+integration(rad-cal_radius,rad+cal_radius,inte_interval));
		}
		else if (cal_radius==0)
		{
			data[i] = G0*den*2*pi*(sqrt(rad*rad+t_dep*t_dep)-sqrt(rad*rad+b_dep*b_dep)+b_dep-t_dep);
		}
	}
	return 0;
}

int vec_cy::datout()
{
	ofstream outfile(outfile_name);
	for (int i = 0; i < M*N; ++i)
	{
		outfile<<xmin+x_interval*(i%N)<<" "<<ymin+y_interval*(i/N)<<" "<<setprecision(16)<<data[i]<<endl;
	}
	outfile.close();
	return 0;
}