// FitCDlg.cpp : implementation file
//
#include "stdafx.h"
#include "FitC.h"
#include "FitCDlg.h"
#include "math.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define pi 3.14159265359
struct FXY{double x,y;};
struct FXYZ{double x,y,z;};
int all_gauss( int n, double *a, double *b);
/////////////////////////////////////////////////////////////////////////////
// CFitCDlg dialog
CFitCDlg::CFitCDlg(CWnd* pParent /*=NULL*/)
: CDialog(CFitCDlg::IDD, pParent)
{
//{{AFX_DATA_INIT(CFitCDlg)
// NOTE: the ClassWizard will add member initialization here
//}}AFX_DATA_INIT
m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}
void CFitCDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
//{{AFX_DATA_MAP(CFitCDlg)
DDX_Control(pDX, IDC_LIST1, m_lst);
//}}AFX_DATA_MAP
}
BEGIN_MESSAGE_MAP(CFitCDlg, CDialog)
//{{AFX_MSG_MAP(CFitCDlg)
ON_WM_PAINT()
ON_WM_QUERYDRAGICON()
ON_BN_CLICKED(IDC_BUTTON1, OnButton1)
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CFitCDlg message handlers
BOOL CFitCDlg::OnInitDialog()
{
CDialog::OnInitDialog();
SetIcon(m_hIcon, TRUE); // Set big icon
SetIcon(m_hIcon, FALSE); // Set small icon
// TODO: Add extra initialization here
return TRUE; // return TRUE unless you set the focus to a control
}
// If you add a minimize button to your dialog, you will need the code below
// to draw the icon. For MFC applications using the document/view model,
// this is automatically done for you by the framework.
void CFitCDlg::OnPaint()
{
if (IsIconic())
{
CPaintDC dc(this); // device context for painting
SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);
// Center icon in client rectangle
int cxIcon = GetSystemMetrics(SM_CXICON);
int cyIcon = GetSystemMetrics(SM_CYICON);
CRect rect;
GetClientRect(&rect);
int x = (rect.Width() - cxIcon + 1) / 2;
int y = (rect.Height() - cyIcon + 1) / 2;
// Draw the icon
dc.DrawIcon(x, y, m_hIcon);
}
else
{
CDialog::OnPaint();
}
}
HCURSOR CFitCDlg::OnQueryDragIcon()
{
return (HCURSOR) m_hIcon;
}
void CFitCDlg::OnButton1()
{
// TODO: Add your control notification handler code here
FXY F3[360];
FXYZ p;
double x0,y0,r0,af;
double x,y,z,XX[12];
double a,b,c;
int k,ret;
char st[80];
//造一套数据
FillMemory(F3,sizeof(FXY)*360,0);
m_lst.AddString("造一套数据");
m_lst.AddString("圆心={100,40};半径=300");
x0=100;y0=40;r0=300;
m_lst.AddString("x0=100;y0=40;r0=300");
m_lst.AddString("");
for(k=0;k<360;k++)
{
af=k*pi/180.;p.z=0.;;
p.x=x=r0*cos(af)+x0;
p.y=y=r0*sin(af)+y0;
F3[k].x=p.x;
F3[k].y=p.y;
sprintf(st,"%4d %9.4f %9.4f",k,x,y);m_lst.AddString(st);
}
m_lst.AddString("");
m_lst.AddString("====================");
m_lst.AddString("计算结果");
//进行计算
FillMemory(XX,sizeof(double)*12,0);
for(k=0;k<360;k++) //状态方程式z=a*x+b*y+c
{
x=F3[k].x;
y=F3[k].y;
z=x*x+y*y;
XX[0]+=x*x;
XX[4]+=y*y;
XX[8]+=1;
XX[1]+=x*y;
XX[2]+=x;
XX[5]+=y;
XX[9]+=x*z;XX[10]+=y*z;XX[11]+=z;
}
XX[3]=XX[1];XX[6]=XX[2];XX[7]=XX[5];
ret=all_gauss(3,XX,&XX[9]);
if(ret==0){delete [] XX;return;}
a=XX[9];b=XX[10];c=XX[11];
sprintf(st,"a=%9.4f b=%9.4f c=%9.4f",a,b,c);m_lst.AddString(st);
x0=a/2;y0=b/2;
r0=sqrt(a*a+b*b+4*c)/2.;
m_lst.AddString("圆心和半径:");
sprintf(st,"x0=%9.4f y0=%9.4f r0=%9.4f",x0,y0,r0);m_lst.AddString(st);
}
int all_gauss( int n, double *a, double *b)
{
int l,k,i,j,is,p,q;//js 用于记忆列交换信息的空间
double d,t;
int *js = (int *) new int[n];
l = 1 ;
for(k=0; k<n-1; k++)
{
d = 0.0;
for(i=k; i<n; i++)for(j=k; j<n; j++)
{
t = fabs(*(a+i*n+j));
if(t>d){ d = t; js[k] = j; is = i;}
} //选主元
if(d+1.0==1.0){return 0;} //奇异标志
else
{
if(js[k]!=k)
for(i=0; i<n; i++)
{ //列交换
p = k;
q = js[k];
t = *(a+i*n+p);
*(a+i*n+p)=*(a+i*n+q);
*(a+i*n+q)=t;
}
if(is!=k)
{
for(j=k; j<n; j++)//行交换
{
t = *(a+k*n+j);
*(a+k*n+j)=*(a+is*n+j);
*(a+is*n+j) = t;
}
t = b[k];
b[k] = b[is];
b[is] = t;
}
}
d = *(a+k*n+k);
for(j=k+1; j<n; j++ )
{
*(a+k*n+j) /= d ;
}
b[k]/=d;
for(i=k+1; i<n; i++)
{ //消元
for(j=k+1; j<n; j++)
{
p = j;
*(a+i*n+p)-= (*(a+i*n+k)) * (*(a+k*n+j));
}
b[i] -= (*(a+i*n+k))*b[k];
}
}
d =*(a+(n-1)*n+n-1);
if(fabs(d)+1.0==1.0){delete [] js;return 0;}
b[n-1]/=d;
for(i=n-2; i>=0; i--)
{
t = 0.0;
for(j=i+1; j<n; j++) t += *(a+i*n+j)*b[j];
b[i] -= t;
}
js[n-1] = n-1;
for(k=n-1; k>=0; k--)if(js[k]!=k)
{
t = b[k];
b[k] = b[js[k]];
b[js[k]] = t; //恢复解向量
}
delete [] js;
return 1;
} // all_gauss