// em.cpp: implementation of the em class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "cp.h"
#include "em.h"
#include "math.h"
#include "stdlib.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

em::em()
{
	pi=3.1415926;
	Ms=1780;
	r=0.0028;
	He=3000;
}

em::~em()
{

}

double em::func0(double st0)
{
	double x0=0;
	double x1=pi/4;
	double f1,f2;
	double sf=0.0001;
	for(;;)
	{
		if(fabs(x0-x1)<sf)break;
		x0=x1;
		f1=(Ms-ha)*sin(2*x0)-2*He*sin(x0-st0);
		f2=(Ms-ha)*2*cos(2*x0)-2*He*cos(x0-st0);
		x1=x0-f1/f2;
	}
	return x1;
}

double em::cpvar(double Ku)
{
	fm=r*Ms;
	ha=Ku*139.6*1800/Ms;
	fa=r*ha;
	return 0;
}

void em::draw(CPaintDC *pdc,Input*pinput)
{
	double f1,f2,f3,f21,f22,f23,u22,u33,u32p;
	int x0=50,y0=450;
	double x,y1,y2,y3,y21,y22,y23,x1;
	char c[30];
	CString cs;
	cpvar(pinput->m_Ku);
	int aa=7,bb=30;
	pdc->MoveTo(x0,20);
	pdc->LineTo(x0,y0);
	pdc->LineTo(x0+700,y0);
	for(int j=0;j<=90;j+=10)
	{
		pdc->MoveTo(x0+j*aa,y0);
		pdc->LineTo(x0+j*aa,y0-10);
		pdc->TextOut(x0+j*aa-6,y0+2,itoa(j,c,10));
	}
	for(j=1;j<=12;j+=1)
	{
		pdc->MoveTo(x0,-30*j+y0);
		pdc->LineTo(x0+10,-30*j+y0);
		pdc->TextOut(x0-20,-30*j+y0-7,itoa(j,c,10));
	}

	double st0=0;
	st0=st0/180*pi;
	f0=getf0(st0);
	f1=sqrt(f0*(f0-fa*pow2(sin(st))));
	f2=sqrt(f0*(f0+fm*pow2(sin(st))-fa*pow2(sin(st))));
	f3=sqrt(f0*(f0+fm-fa*pow2(sin(st))));
	f21=sqrt(f0*f0+f0*(fm-fa)*pow2(sin(st)));
	f22=sqrt((f0+fm)*(f0-fa*pow2(sin(st))));
	if(pow2(tan(st))<(fa/f0))
		f23=sqrt((f0+fm)*(f0-fa*pow2(sin(st))));
	else f23=sqrt(pow2(f0)+f0*(fm-fa)*pow2(sin(st)));
	x=st0*7+x0;
	y1=y0-f1*30;
	y2=y0-f2*30;
	y3=y0-f3*30;
	y21=y0-f21*30;
	y22=y0-f22*30;
	y23=y0-f23*30;


	for(double i=0;i<=90;i+=0.01)
	{
		st0=i;
		st0=st0/180*pi;
		f0=getf0(st0);
		f1=sqrt(f0*(f0-fa*pow2(sin(st))));
		f2=sqrt(f0*(f0+fm*pow2(sin(st))-fa*pow2(sin(st))));
		f3=sqrt(f0*(f0+fm-fa*pow2(sin(st))));
		f21=sqrt(f0*f0+f0*(fm-fa)*pow2(sin(st)));
		f22=sqrt((f0+fm)*(f0-fa*pow2(sin(st))));
		if(f21<f22){f23=f22;}			
		else {f23=f21;}

		x1=x;
		x=i*7+x0;
		if(pinput->m_f1==TRUE)
		{
			y1=y0-f1*30;
			pdc->SetPixel(x,y1,RGB(0,0,0));
		}
		if(pinput->m_f2==TRUE)
		{
			y2=y0-f2*30;
			pdc->SetPixel(x,y2,RGB(0,0,0));
		}
		if(pinput->m_f3==TRUE)
		{
			y3=y0-f3*30;
			pdc->SetPixel(x,y3,RGB(0,0,0));
		}
		if(pinput->m_f4==TRUE&&f22>f21)
		{
			y21=y0-f21*30;
			pdc->SetPixel(x,y21,RGB(22,0,0));
		}
		if(pinput->m_f5==TRUE&&f22>f21)
		{
			y22=y0-f22*30;
			pdc->SetPixel(x,y22,RGB(0,0,0));
		}
		int m=0;
		if(pinput->m_f7==TRUE)
		{
			double t,p;
			f=f23+0.01;
			m=1;
			int i=0;
			for(;;)
			{
				u22=getu22(st);
				u33=getu33(st);
				u32p=getu32p(st);
				t=u22/u33;
				//t=getp2(st);
				if(t<=0){break;}
				p=sqrt(t);
				t=u33*p;
				p=pow2(t)+1+u32p;
				if((p-2*t)/(p+2*t)<=1){break;}
				f=f+0.01;
				i++;
				m=0;
			}
			t=y0-f*30;
			//char c[20];
			//pdc->TextOut(x,t-20,itoa(m,c,10));
			if(i>0)pdc->SetPixel(x,t,RGB(0,0,0));
		}
		if(pinput->m_f6==TRUE)
		{
			y23=y0-f23*30;
			if(m==0)pdc->SetPixel(x,y23,RGB(0,0,0));
		}
	}

	if(pinput->m_f1==TRUE)pdc->TextOut(x+20,y1,"f1");
	if(pinput->m_f2==TRUE)pdc->TextOut(x+20,y2,"f2");
	if(pinput->m_f3==TRUE)pdc->TextOut(x+20,y3,"f3");
	if(pinput->m_f4==TRUE)pdc->TextOut(x+20,y21,"f4");
	if(pinput->m_f5==TRUE)pdc->TextOut(x+20,y22,"f5");
	if(pinput->m_f6==TRUE)pdc->TextOut(x+20,y22,"f6");

	cs=gcvt(st*180/pi,10,c);
	//dc->TextOut(100,160,cs);
	cs=gcvt(st0*180/pi,10,c);
	//dc->TextOut(200,160,cs);
	cs=gcvt(f1,10,c);
	//pdc->TextOut(0,0,cs);
	cs=gcvt(f2,10,c);
	//pdc->TextOut(0,40,cs);
	cs=gcvt(f3,10,c);
	//pdc->TextOut(0,80,cs);
	cs=gcvt(f21,10,c);
	//pdc->TextOut(0,120,cs);
	cs=gcvt(f22,10,c);
	//pdc->TextOut(0,160,cs);
}

double em::getf0(double st0)
{
	if(fabs(st0)<0.001){st=0;return r*(He-Ms+ha);}
	if(fabs(st0-pi/2)<0.001){st=90;return r*He;}
	st=func0(st0);
	return r*He*sin(st0)/sin(st);
}

double em::pow2(double x)
{
	return x*x;
}

double em::getu22(double st)
{
	double t1,t2;
	t1=f0-fa*pow2(sin(st));
	t2=t1*f0-pow2(f);
	t1=t1*fm+t2;
	return t1/t2;
}

double em::getu32p(double st)
{
	return -pow2(f*fm*sin(st)/(f0*(f0-fa*pow2(sin(st)))-pow2(f)));
}

double em::getu33(double st)
{
	double t;
	t=f0*(f0-fa*pow2(sin(st)))-pow2(f);
	return (t+f0*fm*pow2(sin(st)))/t;
}

double em::getp2(double st)
{
	return(((f0+fm)*(f0-fa*pow2(sin(st)))-f*f)/(f0*f0+f0*(fm-fa)*pow2(sin(st)))-f*f);
}