// kawView.cpp : CKawView クラスの動作の定義を行います。
//

#include "stdafx.h"
#include "kaw.h"

#include "kawDoc.h"
#include "kawView.h"

#include "myDialog.h"
#include "myDialog2.h"
#include "myDialog4.h"
#include "myDialog5.h"
#include "myDialog6.h"

#include<math.h>
#include<fstream.h>
#include<conio.h>

#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
#define NR_END 1
#define FREE_ARG char*

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

	void nrerror(char error_text[]);
	double *vector(long nl, long nh);
	int *ivector(long nl, long nh);
	double **matrix(long nrl, long nrh, long ncl, long nch);
	void free_vector(double *v, long nl, long nh);
	void free_ivector(int *v, long nl, long nh);
	void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);


	void mrqmin(double x[], double y[], double sig[], int ndata, double a[], int ma,
				int lista[], int mfit, double **covar, double **alpha, double *chisq,
				void (*funcs)(double, double [], double *, double [], int), double *alamda);

	void fgauss(double x, double a[], double *y, double dyda[], int na);
	void fgauss2(double x, double a[], double *y, int na);

	int selGaus;
	double Lattice, Tube;

	double bessj0(double x);
	double atom(double s);

/////////////////////////////////////////////////////////////////////////////
// CKawView

IMPLEMENT_DYNCREATE(CKawView, CView)

BEGIN_MESSAGE_MAP(CKawView, CView)
	//{{AFX_MSG_MAP(CKawView)
	ON_COMMAND(ID_FILE_OPEN, OnFileOpen)
	ON_COMMAND(ID_FILE_SAVE_AS, OnFileSaveAs)
	ON_COMMAND(ID_FILE_SAVE, OnFileSave)
	ON_COMMAND(ID_Region, OnRegion)
	ON_COMMAND(ID_Peak, OnPeak)
	ON_COMMAND(ID_ResetRegion, OnResetRegion)
	ON_COMMAND(ID_Marq, OnMarq)
	ON_COMMAND(ID_SaveMarq, OnSaveMarq)
	ON_WM_LBUTTONDOWN()
	ON_WM_LBUTTONUP()
	ON_WM_MOUSEMOVE()
	ON_WM_RBUTTONDOWN()
	ON_WM_LBUTTONDBLCLK()
	ON_UPDATE_COMMAND_UI(ID_Region, OnUpdateRegion)
	ON_UPDATE_COMMAND_UI(ID_Peak, OnUpdatePeak)
	ON_UPDATE_COMMAND_UI(ID_SaveMarq, OnUpdateSaveMarq)
	ON_UPDATE_COMMAND_UI(ID_Marq, OnUpdateMarq)
	ON_WM_KEYDOWN()
	ON_COMMAND(ID_LS_Condition, OnLSCondition)
	ON_COMMAND(ID_RegionInput, OnRegionInput)
	ON_COMMAND(ID_CalPattern, OnCalPattern)
	ON_COMMAND(ID_seqMarq, OnseqMarq)
	ON_UPDATE_COMMAND_UI(ID_seqMarq, OnUpdateseqMarq)
	ON_COMMAND(ID_bg, Onbg)
	ON_UPDATE_COMMAND_UI(ID_bg, OnUpdatebg)
	//}}AFX_MSG_MAP
	// 標準印刷コマンド
	ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CKawView クラスの構築/消滅

CKawView::CKawView()
{
	// TODO: この場所に構築用のコードを追加してください。
	Marq=0;
	SetPeak=0;
	MRegion=2;
	Xmin=0;
	Xmax=0;
	Ymin=0;
	Ymax=0;
	NumPeak=0;CPeak=0;
	Wcoef[0]=1.0;Wcoef[1]=1.0;Wcoef[2]=1.0;Wcoef[3]=1.0;Wcoef[4]=1.0;
	selGaus=1;
	Convergence=0.001;
    Xcolumn=1; Ycolumn=2; SkipLine=0; Lastline=9999;
	Wave=0.61821;
	Lattice=17.4, Tube=7.1, FWHM=0.002, FWHM_v=0.0, BackA=0.0, BackB=0.0, BackC=0.0, Scale=10.0, Thermal=10.0;
	F_calc=0;
	Lattice_min=15.0, Lattice_max=18.0, Lattice_step=0.1;
	Tube_min=6.0, Tube_max=8.0, Tube_step=0.1;
	F_seqMarq=0;
}

CKawView::~CKawView()
{
}

BOOL CKawView::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: この位置で CREATESTRUCT cs を修正して Window クラスまたはスタイルを
	//  修正してください。

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CKawView クラスの描画

void CKawView::OnDraw(CDC* pDC)
{	
	double arg;
	double Yc;

	CKawDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	// TODO: この場所にネイティブ データ用の描画コードを追加します。

	int xt,yt, i, j;
	int firstY=0;
	double Yi,Xi;

    pDC->TextOut(10,10, openfilename);

    for (i=1; i<=pDoc->ndata; i++){

		Xi=pDoc->x[i];
		Yi=pDoc->y[i];
		if (Xi>=Xmin && Xi<=Xmax &&
			Yi>=Ymin && Yi<=Ymax) {
			xt=(int) 550*(Xi-Xmin)/(Xmax-Xmin)+100;
			yt=(int) -300*(Yi-Ymin)/(Ymax-Ymin)+350;

			pDC->SetPixel(xt,yt,RGB(0,0,255));
			pDC->SetPixel(xt+1,yt,RGB(0,0,255));
			pDC->SetPixel(xt,yt+1,RGB(0,0,255));
			pDC->SetPixel(xt+1,yt+1,RGB(0,0,255));
		}

		if (F_calc==1) {
 			if (Xcalc[i]>=Xmin && Xcalc[i]<=Xmax &&
				Ycalc[i]>=Ymin && Ycalc[i]<=Ymax) {
				xt=(int) 550*(Xcalc[i]-Xmin)/(Xmax-Xmin)+100;
				yt=(int) -300*(Ycalc[i]-Ymin)/(Ymax-Ymin)+350;

				pDC->SetPixel(xt,yt,RGB(255,0,255));
				pDC->SetPixel(xt+1,yt,RGB(255,0,255));
				pDC->SetPixel(xt,yt+1,RGB(255,0,255));
				pDC->SetPixel(xt+1,yt+1,RGB(255,0,255));
			}
		}

	}

    char name[20], *t1;

    if (pDoc->ndata>0) {
        pDC->TextOut(250,30, "Curve Fitting");

        t1 = gcvt(Xmin, 8, name);
        pDC->TextOut(100,360, t1);
        t1 = gcvt(Xmax, 8, name);
        pDC->TextOut(640,360, t1);

        pDC->SetTextColor(RGB(255,0,0));
        t1 = gcvt(Ymax, 8, name);
        pDC->TextOut(35,50, t1);
        t1 = gcvt(Ymin, 8, name);
        pDC->TextOut(35,330, t1);
    }


	if (pDoc->ndata>0) {
		CBrush NewBrush;
		NewBrush.CreateSolidBrush(RGB(125,125,125));
		CBrush* OldBrush=pDC->SelectObject(&NewBrush);

		CRect DRect(100,50,650,350);
		pDC->FrameRect(DRect, &NewBrush);

		pDC->SelectObject(OldBrush);
		NewBrush.DeleteObject();
	}
    
	if (F_seqMarq==1) {
	   pDC->TextOut(660,30, "Now calculating");
	}

}

////////////////////////////////////
double bessj0(double x)
{
	double ax,z;
	double xx,y,ans,ans1,ans2;

	if ((ax=fabs(x)) < 8.0) {
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
			+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
			+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans=ans1/ans2;
	} else {
		z=8.0/ax;
		y=z*z;
		xx=ax-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}

double atom(double s)
{
	double a[5],b[5],c;
	double ans;

	a[1]=2.31; a[2]=1.02; a[3]=1.5886; a[4]=0.865;
	b[1]=20.8439; b[2]=10.2075; b[3]=0.5687; b[4]=51.6512;
	c=0.2156;

	ans=0.0;

	for (int i=1; i<=4; i++){
		ans=ans+a[i]*exp(-b[i]*s*s);
	}

	ans=ans+c;

	return ans;
}


/////////////////////////////////////////////////////////////////////////////
// CKawView クラスの印刷

BOOL CKawView::OnPreparePrinting(CPrintInfo* pInfo)
{
	// デフォルトの印刷準備
	return DoPreparePrinting(pInfo);
}

void CKawView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: 印刷前の特別な初期化処理を追加してください。
}

void CKawView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: 印刷後の後処理を追加してください。
}

/////////////////////////////////////////////////////////////////////////////
// CKawView クラスの診断

#ifdef _DEBUG
void CKawView::AssertValid() const
{
	CView::AssertValid();
}

void CKawView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CKawDoc* CKawView::GetDocument() // 非デバッグ バージョンはインラインです。
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CKawDoc)));
	return (CKawDoc*)m_pDocument;
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CKawView クラスのメッセージ ハンドラ

void CKawView::OnFileOpen() 
{
	
	CKawDoc* pDoc = GetDocument();

	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	CFileDialog filedlg(TRUE, "", "*.*");
	CString filename;
	FILE *ifp;
	char s[256];
	int n,m,k,nd;

	if (filedlg.DoModal() == IDOK) {

		filename = filedlg.GetPathName();
        openfilename=filename;
        ifp = fopen(filename, "r");
	}

	CMyDialog4 mD;

	int ret4;

	mD.m_x_column=Xcolumn;
	mD.m_y_column=Ycolumn;
	mD.m_skipline=SkipLine;
	mD.m_lastline=Lastline;
	mD.m_wavelength=Wave;

    ret4=mD.DoModal();

	if (ret4==IDOK) {
	    Xcolumn=mD.m_x_column;
	    Ycolumn=mD.m_y_column;
		SkipLine=mD.m_skipline;
		Lastline=mD.m_lastline;
	    Wave=mD.m_wavelength;
	}

    for (n=0; n<SkipLine; n++) {
         fgets(s,256,ifp);
	}
        
    nd=0;
        
    while (fgets(s,256,ifp)!=NULL){

        m=0;
        k=0;

		if (nd>=Lastline) break;
         
        if (strlen(s)<=1) break;

        if (strlen(s)) {

            for(n=0; n<strlen(s); n++) {

                if (k>=2) {
                    m=0;
                    k=0;

                    break;
                }
       
                if (isspace(s[n]) || s[n]==',') {
                    m=0;
                }
                else {
                    if (m==0) {
                       if (k==0) { nd=nd+1;}
                        k=k+1;
                       if (k==Xcolumn) {
                           pDoc->x[nd]=atof(&s[n]);
                       }
                       else if (k==Ycolumn) {
                           pDoc->y[nd]=atof(&s[n]);
                       }
                    }
                    m=1;    
               }
            } 
       }

	}

	fclose(ifp);
       
	pDoc->ndata=nd;

    for (int i=1; i<=nd; i++) {
         pDoc->x[i]=sin(pDoc->x[i]/2.0*3.14159/180.0)/Wave;
	}

    pDoc->startA=pDoc->x[1];
    pDoc->stepA=pDoc->x[2]-pDoc->x[1];
    pDoc->stopA=pDoc->x[pDoc->ndata];

	pDoc->dMax=pDoc->y[1];
	pDoc->dMin=pDoc->y[1];

    for (i=1; i<=pDoc->ndata; i++) {

		if (pDoc->dMax < pDoc->y[i]) {
			pDoc->dMax = pDoc->y[i];
		}

		if (pDoc->dMin > pDoc->y[i]) {
			pDoc->dMin = pDoc->y[i];
		}
	}

    Xmin=pDoc->x[1];
	Xmax=pDoc->x[pDoc->ndata];

	if (Xmin>Xmax) {
        Xmax=pDoc->x[1];
	    Xmin=pDoc->x[pDoc->ndata];
	}
	
	Ymin=pDoc->dMin;
	Ymax=pDoc->dMax;

    Marq=0;
	MRegion=2;

    RedrawWindow();	

}

void CKawView::OnFileSaveAs() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	CKawDoc* pDoc = GetDocument();

	CFileDialog filedlg(FALSE, "dat", "*.dat");
	CString filename;
	FILE *ofp;
	int m=0,nd=0;
	
	if (filedlg.DoModal() == IDOK) {

		filename = filedlg.GetPathName();

		ofp = fopen(filename, "w");

		for (int i=0; i<4; i++) {
	        fprintf(ofp, "\n","");
		}

		fprintf(ofp, "%8.3f%8.3f%8.3f%10s\n",pDoc->startA,pDoc->stepA,pDoc->stopA,"from RINT");

		while (nd<pDoc->ndata){

			m=m+1;

			if (m==11) {
				fprintf(ofp, "\n","");
				m=0;
			}
			else {
				nd=nd+1;
				fprintf(ofp, "%8d", (int)pDoc->y[nd]);
			}
		}

		fprintf(ofp, "\n","");

		fclose(ofp);
	}

}

void CKawView::OnFileSave() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	CKawDoc* pDoc = GetDocument();

	CFileDialog filedlg(FALSE, "dat", "*.dat");
	CString filename;
	FILE *ofp;
	int m=0,nd=0;
	
	if (filedlg.DoModal() == IDOK) {

		filename = filedlg.GetPathName();

		ofp = fopen(filename, "w");

//		fprintf(ofp, "%d\n", pDoc->ndata);

		for (int i=1; i<=pDoc->ndata; i++){
			fprintf(ofp, "%f %f\n", pDoc->x[i], pDoc->y[i]);
		}

		fclose(ofp);
	}

}


void mrqmin(double x[], double y[], double sig[], int ndata, double a[], int ma,
	int lista[], int mfit, double **covar, double **alpha, double *chisq,
	void (*funcs)(double, double [], double *, double [], int), double *alamda)
{

	void covsrt(double **covar, int ma, int lista[], int mfit);

	void gaussj(double **a, int n, double **b, int m);

	void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
		int ma, int lista[], int mfit, double **alpha, double beta[],
		double *chisq, void (*funcs)(double,double *,double *,double *,int));


	int k,kk,j,ihit;
	static double *da,*atry,**oneda,*beta,ochisq;

	if (*alamda < 0.0) {
		oneda=matrix(1,mfit,1,1);
		atry=vector(1,ma);
		da=vector(1,ma);
		beta=vector(1,ma);
		kk=mfit+1;
		for (j=1;j<=ma;j++) {
			ihit=0;
			for (k=1;k<=mfit;k++)
				if (lista[k] == j) ihit++;
			if (ihit == 0)
				lista[kk++]=j;
			else if (ihit > 1) nrerror("Bad lista permutation in mrqmin-1");
		}
		if (kk != ma+1) nrerror("Bad lista permutation in mrqmin-2");
		*alamda=0.001;
		mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,chisq,funcs);
		ochisq=(*chisq);
	}
	for (j=1;j<=mfit;j++) {
		for (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k];
		covar[j][j]=alpha[j][j]*(1.0+(*alamda));
		oneda[j][1]=beta[j];
	}
	gaussj(covar,mfit,oneda,1);
	for (j=1;j<=mfit;j++)
		da[j]=oneda[j][1];
	if (*alamda == 0.0) {
		covsrt(covar,ma,lista,mfit);
		free_vector(beta,1,ma);
		free_vector(da,1,ma);
		free_vector(atry,1,ma);
		free_matrix(oneda,1,mfit,1,1);
		return;
	}
	for (j=1;j<=ma;j++) atry[j]=a[j];
	for (j=1;j<=mfit;j++)
		atry[lista[j]] = a[lista[j]]+da[j];
	mrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,chisq,funcs);
	if (*chisq < ochisq) {
		*alamda *= 0.1;
		ochisq=(*chisq);
		for (j=1;j<=mfit;j++) {
			for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k];
			beta[j]=da[j];
			a[lista[j]]=atry[lista[j]];
		}
	} else {
		*alamda *= 10.0;
		*chisq=ochisq;
	}
	return;
}

void covsrt(double **covar, int ma, int lista[], int mfit)
{
	int i,j;
	double swap;

	for (j=1;j<ma;j++)
		for (i=j+1;i<=ma;i++) covar[i][j]=0.0;
	for (i=1;i<mfit;i++)
		for (j=i+1;j<=mfit;j++) {
			if (lista[j] > lista[i])
				covar[lista[j]][lista[i]]=covar[i][j];
			else
				covar[lista[i]][lista[j]]=covar[i][j];
		}
	swap=covar[1][1];
	for (j=1;j<=ma;j++) {
		covar[1][j]=covar[j][j];
		covar[j][j]=0.0;
	}
	covar[lista[1]][lista[1]]=swap;
	for (j=2;j<=mfit;j++) covar[lista[j]][lista[j]]=covar[1][j];
	for (j=2;j<=ma;j++)
		for (i=1;i<=j-1;i++) covar[i][j]=covar[j][i];
}


void gaussj(double **a, int n, double **b, int m)
{
	int *indxc,*indxr,*ipiv;
	int i,icol,irow,j,k,l,ll;
	double big,dum,pivinv,temp;

	indxc=ivector(1,n);
	indxr=ivector(1,n);
	ipiv=ivector(1,n);
	for (j=1;j<=n;j++) ipiv[j]=0;
	for (i=1;i<=n;i++) {
		big=0.0;
		for (j=1;j<=n;j++)
			if (ipiv[j] != 1)
				for (k=1;k<=n;k++) {
					if (ipiv[k] == 0) {
						if (fabs(a[j][k]) >= big) {
							big=fabs(a[j][k]);
							irow=j;
							icol=k;
						}
					} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
				}
		++(ipiv[icol]);
		if (irow != icol) {
			for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
			for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
		}
		indxr[i]=irow;
		indxc[i]=icol;
		if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
		pivinv=1.0/a[icol][icol];
		a[icol][icol]=1.0;
		for (l=1;l<=n;l++) a[icol][l] *= pivinv;
		for (l=1;l<=m;l++) b[icol][l] *= pivinv;
		for (ll=1;ll<=n;ll++)
			if (ll != icol) {
				dum=a[ll][icol];
				a[ll][icol]=0.0;
				for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
				for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
			}
	}
	for (l=n;l>=1;l--) {
		if (indxr[l] != indxc[l])
			for (k=1;k<=n;k++)
				SWAP(a[k][indxr[l]],a[k][indxc[l]]);
	}
	free_ivector(ipiv,1,n);
	free_ivector(indxr,1,n);
	free_ivector(indxc,1,n);
}

void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
	int ma, int lista[], int mfit, double **alpha, double beta[],
	double *chisq, void (*funcs)(double,double *,double *,double *,int))
{
	int k,j,i;
	double ymod,wt,sig2i,dy,*dyda;

	dyda=vector(1,ma);
	for (j=1;j<=mfit;j++) {
		for (k=1;k<=j;k++) alpha[j][k]=0.0;
		beta[j]=0.0;
	}
	*chisq=0.0;
	for (i=1;i<=ndata;i++) {
		(*funcs)(x[i],a,&ymod,dyda,ma);
		sig2i=1.0/(sig[i]*sig[i]);
		dy=y[i]-ymod;
		for (j=1;j<=mfit;j++) {
			wt=dyda[lista[j]]*sig2i;
			for (k=1;k<=j;k++)
				alpha[j][k] += wt*dyda[lista[k]];
			beta[j] += dy*wt;
		}
		(*chisq) += dy*dy*sig2i;
	}
	for (j=2;j<=mfit;j++)
		for (k=1;k<=j-1;k++) alpha[k][j]=alpha[j][k];
	free_vector(dyda,1,ma);
}

void fgauss(double x, double a[], double *y, double dyda[], int na)
{
    double d_Value[5],Multi[5];

    d_Value[1]=sqrt(3)*0.5*Lattice;
	d_Value[2]=sqrt(3)*0.5*Lattice/sqrt(3);
	d_Value[3]=sqrt(3)*0.5*Lattice/sqrt(4);
	d_Value[4]=sqrt(3)*0.5*Lattice/sqrt(7);
	Multi[1]=1.0;
	Multi[2]=1.0;
	Multi[3]=1.0;
	Multi[4]=2.0;

	double arg, Scatter;
	double y1, y2, y3;
	double W;

	y1=0.0;
	y2=0.0;
	y3=0.0;
	*y=0.0;

	for (int j=1; j<=4; j++) {


            W=a[2]+a[3]*x;
			arg=(x-1.0/(2.0*d_Value[j]))/W;
			y1=y1+Multi[j]*exp(-arg*arg);
			y2=y2+Multi[j]*2.0*(x-1.0/(2.0*d_Value[j]))*(x-1.0/(2.0*d_Value[j]))/(W*W*W)*exp(-arg*arg);
			y3=y3+Multi[j]*2.0*(x-1.0/(2.0*d_Value[j]))*(x-1.0/(2.0*d_Value[j]))*x/(W*W*W)*exp(-arg*arg);

	}

	Scatter=atom(x)*bessj0(4.0*3.14159*x*Tube);

	*y += a[1]*y1*Scatter*Scatter*exp(-a[4]*x*x)+a[5]*x*x+a[6]*x+a[7];

	dyda[1]=y1*Scatter*Scatter*exp(-a[4]*x*x);
	dyda[2]=a[1]*y2*Scatter*Scatter*exp(-a[4]*x*x);
	dyda[3]=a[1]*y3*Scatter*Scatter*exp(-a[4]*x*x);
	dyda[4]=a[1]*y1*Scatter*Scatter*(-x*x)*exp(-a[4]*x*x);
	dyda[5]=x*x;
	dyda[6]=x;
	dyda[7]=1.0;

}


void fgauss2(double x, double a[], double *y, int na)
{
	int i;
	double ex,arg;

	*y=0.0;
	if (selGaus==0){
		for (i=1;i<=na-1;i+=3) {
			arg=(x-a[i+1])/a[i+2];
			ex=exp(-arg*arg);
			*y += a[i]*ex;
		}
	}
	else  if (selGaus==1) {
		for (i=1;i<=na-1;i+=3) {
			arg=(x-a[i+1])/a[i+2];
			ex=1.0/(arg*arg+1.0);
			*y += a[i]*ex;
		}
	}

	

}

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
	exit(1);
}

double *vector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl+NR_END;
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;

	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl+NR_END;
}

double **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

void free_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by matrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}



void CKawView::OnRegion() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	SetPeak=0;
	MRegion=0;
	NumPeak=0;
	Marq=0;
}

void CKawView::OnPeak() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	SetPeak=1;
}

void CKawView::OnResetRegion() 
{
	CKawDoc* pDoc = GetDocument();

	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	SetPeak=0;
	MRegion=0;
	Marq=0;
	NumPeak=0;

    Xmin=pDoc->x[1];
	Xmax=pDoc->x[pDoc->ndata];
	Ymin=pDoc->dMin;
	Ymax=pDoc->dMax;

	RedrawWindow();	

}

void CKawView::OnMarq() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	
	CKawDoc* pDoc = GetDocument();

    int i;

    static double x[4096], y[4096];
    static double sig[4096];
    int ndata;
    double a[20];
    int ma;
    int lista[20];
    int mfit;
    double **covar;
    double **alpha;
    double *chisq;
    double *alamda;
    double aa,cc;
	double *chiold, *chinew;
	double Chiinitial, Chiold, Chinew;

	int n;

	n=0;

    for (i=1; i<=pDoc->ndata; i++)
    {

		if (pDoc->x[i]>=Xmin && pDoc->x[i]<=Xmax) {
			n=n+1;
            x[n] = pDoc->x[i];
            y[n] = pDoc->y[i];
            sig[n] = 1.0/x[n];
		}
    }

    a[1]=Scale;
	a[2]=FWHM;
	a[3]=FWHM_v;
	a[4]=Thermal;
	a[5]=BackA;
	a[6]=BackB;
	a[7]=BackC;

    ma=7;
    mfit=7;

	for (i=1; i<=ma; i++){
		lista[i]=i;
	}

    aa=-1.0;cc=0.0;

    alamda=&aa;
    chisq=&cc;

    ndata=n;

    covar=matrix(1,ma,1,ma);
    alpha=matrix(1,ma,1,ma);

    mrqmin(x, y, sig, ndata, a, ma, lista, mfit, 
               covar, alpha, chisq, fgauss, alamda);

	chiold=chisq;
	Chiold=*chiold;

    for (i=1;i<=10;i++)
    {

        mrqmin(x, y, sig, ndata, a, ma, lista, mfit, 
               covar, alpha, chisq, fgauss, alamda);

		chinew=chisq;
		Chinew=*chinew;


		if (sqrt((Chiold-Chinew)*(Chiold-Chinew))<Convergence*(Chiold)) {
			MessageBox("収束しました");
			break;
		}
		else {
			MessageBox("収束していません");
			Chiold=Chinew;
		}

    }


    Marq=1;

    char name[20], *t1;
    
    t1 = gcvt(a[1], 8, name);
    AfxMessageBox(name, MB_OK|MB_ICONEXCLAMATION);

    Scale=a[1];
	FWHM=a[2];
	FWHM_v=a[3];
	Thermal=a[4];
	BackA=a[5];
	BackB=a[6];
	BackC=a[7];

	double Xi, y1, arg, Scatter;
	double d_Value[5], Multi[5];
	double W;

    d_Value[1]=sqrt(3)*0.5*Lattice;
	d_Value[2]=sqrt(3)*0.5*Lattice/sqrt(3);
	d_Value[3]=sqrt(3)*0.5*Lattice/sqrt(4);
	d_Value[4]=sqrt(3)*0.5*Lattice/sqrt(7);
	Multi[1]=1.0;
	Multi[2]=1.0;
	Multi[3]=1.0;
	Multi[4]=2.0;

    for (i=1; i<=pDoc->ndata; i++){

		Xi=pDoc->x[i];

		if (Xi>=Xmin && Xi<=Xmax)

		y1=0.0;

		for (int j=1; j<=4; j++) {

			W=FWHM+FWHM_v*Xi;
			arg=(Xi-1.0/(2.0*d_Value[j]))/W;
			y1=y1+Multi[j]*exp(-arg*arg);
		}

		Scatter=atom(Xi)*bessj0(4.0*3.14159*Xi*Tube);
		Ycalc[i]=Scale*y1*Scatter*Scatter*exp(-Thermal*Xi*Xi)+BackA*Xi*Xi+BackB*Xi+BackC;

	}

    RedrawWindow();	
	
}

void CKawView::OnSaveMarq() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	
}

void CKawView::OnLButtonDown(UINT nFlags, CPoint point) 
{
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
	if (SetPeak==0 && MRegion==0) {
		MRegion=1;
		RegionSX=point.x;
		RegionSY=point.y;
	}
}

void CKawView::OnLButtonUp(UINT nFlags, CPoint point) 
{
	CDC* pDC=GetDC();

	double temp1,temp2;
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
	if (SetPeak==0 && MRegion==1) {
		MRegion=0;
	   
		temp1=(Xmax-Xmin)*(RegionSX-100)/550.0+Xmin;
		temp2=(Xmax-Xmin)*(point.x-100)/550.0+Xmin;
	   
		if (temp1>temp2) {
	   	   Xmin=temp2;
		   Xmax=temp1;
		}
		else {
	   	   Xmin=temp1;
		   Xmax=temp2;
		}

		temp1=(Ymin-Ymax)*(RegionSY-50)/300+Ymax;
		temp2=(Ymin-Ymax)*(point.y-50)/300+Ymax;
		if (temp1>temp2) {
	   	   Ymin=temp2;
		   Ymax=temp1;
		}
		else {
	   	   Ymin=temp1;
		   Ymax=temp2;
		}

	}


    RedrawWindow();	

}

void CKawView::OnMouseMove(UINT nFlags, CPoint point) 
{
	CDC* pDC=GetDC();

	int EX,EY;
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
	if (SetPeak==0 && MRegion==1) {
		EX=point.x;
		EY=point.y;
		RedrawWindow();
		pDC->MoveTo(RegionSX,RegionSY);
		pDC->LineTo(RegionSX,EY);
		pDC->LineTo(EX,EY);
		pDC->LineTo(EX,RegionSY);
		pDC->LineTo(RegionSX,RegionSY);
	}
  }

void CKawView::OnRButtonDown(UINT nFlags, CPoint point) 
{
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
	if (NumPeak>=1) {
		NumPeak=NumPeak-1;
		CPeak=NumPeak;
	}
    RedrawWindow();	

}

void CKawView::OnLButtonDblClk(UINT nFlags, CPoint point) 
{
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
	if (SetPeak==1 && NumPeak<=5) {
		NumPeak=NumPeak+1;
		CPeak=NumPeak;
		PeakP[NumPeak]=(Xmax-Xmin)*(point.x-100)/550.0+Xmin;
        PeakH[NumPeak]=(Ymin-Ymax)*(point.y-50)/300.0+Ymax;
	}
}

void CKawView::OnUpdateRegion(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
	pCmdUI->SetCheck((SetPeak==0 && MRegion<2) ? 1:0);
}

void CKawView::OnUpdatePeak(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
	pCmdUI->SetCheck((SetPeak==1) ? 1:0);	
}

void CKawView::OnUpdateSaveMarq(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
	pCmdUI->Enable((Marq==1) ? TRUE: FALSE);		
}

void CKawView::OnUpdateMarq(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
	pCmdUI->Enable((NumPeak>0) ? TRUE: FALSE);	
}

void CKawView::OnKeyDown(UINT nChar, UINT nRepCnt, UINT nFlags) 
{
	// TODO: この位置にメッセージ ハンドラ用のコードを追加するかまたはデフォルトの処理を呼び出してください
    CClientDC myDC(this);

    if (nChar=='W') {
		if (NumPeak>0) {
			Wcoef[CPeak]=Wcoef[CPeak]*1.05;
		}
    }
    if (nChar=='N') {
		if (NumPeak>0) {
		    Wcoef[CPeak]=Wcoef[CPeak]*0.95;
		}
    }

	if (nChar==VK_RIGHT) {
		if (NumPeak>0) {
			if ((CPeak+1)<=NumPeak){
		        CPeak=CPeak+1;
			}
		}
	}
	if (nChar==VK_LEFT) {
		if (NumPeak>0) {
			if ((CPeak-1)>=1){
		        CPeak=CPeak-1;
			}
		}
	}

    RedrawWindow();	
	
	CView::OnKeyDown(nChar, nRepCnt, nFlags);
}

void CKawView::OnLSCondition() 
{
    int ret;
	static double cv=0.001;

    CmyDialog myDL;

	myDL.m_edit1=cv;

    ret=myDL.DoModal();
    if (ret==IDOK) {
            cv=myDL.m_edit1;
			Convergence=cv;
    }

}


void CKawView::OnRegionInput() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
    int ret;

    CmyDialog2 myDL;

	myDL.m_xmin=Xmin;
	myDL.m_xmax=Xmax;
	myDL.m_ymin=Ymin;
	myDL.m_ymax=Ymax;

    ret=myDL.DoModal();
    if (ret==IDOK) {
            Xmin=myDL.m_xmin;
		    Xmax=myDL.m_xmax;
			Ymin=myDL.m_ymin;
			Ymax=myDL.m_ymax;
    }

    RedrawWindow();	
	
}




void CKawView::OnCalPattern() 
{
	CKawDoc* pDoc = GetDocument();
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください

	CMyDialog5 mD;

	int ret5;

	mD.m_LatticeC=Lattice;
	mD.m_TubeR=Tube;
	mD.m_FWHM=FWHM;
	mD.m_FWHM_v=FWHM_v;
	mD.m_BackA=BackA;
	mD.m_BackB=BackB;
	mD.m_BackC=BackC;
	mD.m_Scale=Scale;
	mD.m_thermal=Thermal;

    ret5=mD.DoModal();

	if (ret5==IDOK) {
		Lattice=mD.m_LatticeC;
		Tube=mD.m_TubeR;
		FWHM=mD.m_FWHM;
		FWHM_v=mD.m_FWHM_v;
		BackA=mD.m_BackA;
		BackB=mD.m_BackB;
		BackC=mD.m_BackC;
		Scale=mD.m_Scale;
		Thermal=mD.m_thermal;
	}

    double d_Value[5], Multi[5];

    d_Value[1]=sqrt(3)*0.5*Lattice;
	d_Value[2]=sqrt(3)*0.5*Lattice/sqrt(3);
	d_Value[3]=sqrt(3)*0.5*Lattice/sqrt(4);
	d_Value[4]=sqrt(3)*0.5*Lattice/sqrt(7);
	Multi[1]=1.0;
	Multi[2]=1.0;
	Multi[3]=1.0;
	Multi[4]=2.0;

	double Xi, arg;
	double Scatter[2000],Scatter2[2000],AtomS[2000];
	double y1;

    for (int i=1; i<=pDoc->ndata; i++){

		Xi=pDoc->x[i];

		if (Xi>=Xmin && Xi<=Xmax) 

		y1=0.0;

		for (int j=1; j<=4; j++) {

			arg=(Xi-1.0/(2.0*d_Value[j]))/(FWHM + FWHM_v*Xi);
			y1=y1+Multi[j]*exp(-arg*arg);
		}

		Scatter[i]=atom(Xi)*bessj0(4.0*3.14159*Xi*Tube);
		AtomS[i]=atom(Xi);
		Scatter2[i]=Scatter[i]*Scatter[i];

		Ycalc[i]=Scale*y1*Scatter[i]*Scatter[i]*exp(-Thermal*Xi*Xi)+BackA*Xi*Xi+BackB*Xi+BackC;
		Xcalc[i]=Xi;
		

	}

	F_calc=1;
	NumPeak=1;

    RedrawWindow();	
/*
	CFileDialog filedlg(FALSE, "dat", "*.dat");

	CString filename;

	if( filedlg.DoModal() == IDOK)
	{
		filename = filedlg.GetPathName();
                ofstream mystream(filename);

                for (int i=1; i<=pDoc->ndata; i++){
                    mystream << Xcalc[i] << ' ' << Ycalc[i] << ' ' << Scatter2[i] << ' ' << AtomS[i] << endl;
                }

		mystream.close();
	}
*/

}

void CKawView::OnseqMarq() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	CKawDoc* pDoc = GetDocument();

	CMyDialog6 mD;

	int ret6;
	int num_L_step, num_R_step;
	int num_L, num_R;

	mD.m_L_min=Lattice_min;
	mD.m_L_max=Lattice_max;
	mD.m_L_step=Lattice_step;
	mD.m_R_min=Tube_min;
	mD.m_R_max=Tube_max;
	mD.m_R_step=Tube_step;

    ret6=mD.DoModal();


	if (ret6==IDOK) {

		Lattice_min=mD.m_L_min;
		Lattice_max=mD.m_L_max;
		Lattice_step=mD.m_L_step;
		Tube_min=mD.m_R_min;
		Tube_max=mD.m_R_max;
		Tube_step=mD.m_R_step;

	}

	num_L_step= (Lattice_max-Lattice_min)/Lattice_step;
	num_R_step= (Tube_max-Tube_min)/Tube_step;

    int i;

    static double x[4096], y[4096];
    static double sig[4096];
    int ndata;
    double a[20];
    int ma;
    int lista[20];
    int mfit;
    double **covar;
    double **alpha;
    double *chisq;
    double *alamda;
    double aa,cc;
	double *chiold, *chinew;
	double Chiinitial, Chiold, Chinew;
	double Chi_min;
	int num_calc, num_min;

	int n;
	int n_marq;

	CFileDialog filedlg(FALSE, "dat", "*.dat");

	CString filename;

	if( filedlg.DoModal() == IDOK)
	{
		filename = filedlg.GetPathName();
                ofstream mystream(filename);


	n=0;

    for (i=1; i<=pDoc->ndata; i++)
    {

		if (pDoc->x[i]>=Xmin && pDoc->x[i]<=Xmax) {
			n=n+1;
            x[n] = pDoc->x[i];
            y[n] = pDoc->y[i];
            sig[n] = 1.0/x[n];
		}
    }

    ndata=n;

    ma=7;
    mfit=7;

	for (i=1; i<=ma; i++){
		lista[i]=i;
	}

	F_seqMarq=1;
    RedrawWindow();	

	Chi_min=0.0;
	num_min=0;
	num_calc=0;

    for (num_L=0; num_L<=num_L_step; num_L++) {

		Lattice=Lattice_min + num_L * Lattice_step;


		if (kbhit()){
			switch (getch()){
				case 0x41: break;
				case 0x42: break;

			}
		}

		for (num_R=0; num_R<=num_R_step; num_R++) {

     		Tube=Tube_min + num_R * Tube_step;

			a[1]=Scale;
			a[2]=FWHM;
			a[3]=FWHM_v;
			a[4]=Thermal;
			a[5]=BackA;
			a[6]=BackB;
			a[7]=BackC;


			aa=-1.0;cc=0.0;

			alamda=&aa;
			chisq=&cc;

			covar=matrix(1,ma,1,ma);
			alpha=matrix(1,ma,1,ma);

			mrqmin(x, y, sig, ndata, a, ma, lista, mfit, 
               covar, alpha, chisq, fgauss, alamda);


			chiold=chisq;
			Chiold=*chiold;
			Chiinitial=Chiold;


			for (i=1;i<=10;i++) {

				mrqmin(x, y, sig, ndata, a, ma, lista, mfit, 
					covar, alpha, chisq, fgauss, alamda);

				chinew=chisq;
				Chinew=*chinew;


				if (sqrt((Chiold-Chinew)*(Chiold-Chinew))<Convergence*(Chiold)) {
					break;
				}
				else {
					Chiold=Chinew;
				}

			}

			num_calc=num_calc+1;

			if (Chi_min==0.0) {
				Chi_min=Chinew;
				num_min=num_calc;
			}
			else if (Chi_min>Chinew) {
				Chi_min=Chinew;
				num_min=num_calc;
			}

			n_marq=i;

			mystream << Lattice << ' ' << Tube << ' ' <<' ' << Chinew << ' ' 
				     << n_marq << ' ' << Chiinitial << ' ' << Chiold << ' ' 
					 << ' ' << a[1] <<' ' << a[2] <<' ' << a[3] <<' ' << a[4] <<' ' 
					 << a[5] <<' ' << a[6] <<' ' << a[7] <<endl;

		}
	}

	mystream << endl;
	mystream << num_min << endl;

	mystream.close();
	
	}

	F_seqMarq=0;
    RedrawWindow();
	
	MessageBox("計算終了");


}

void CKawView::OnUpdateseqMarq(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
		pCmdUI->Enable((NumPeak>0) ? TRUE: FALSE);	
}

void CKawView::Onbg() 
{
	// TODO: この位置にコマンド ハンドラ用のコードを追加してください
	CKawDoc* pDoc = GetDocument();
	CFileDialog filedlg(FALSE, "dat", "*.dat");

	CString filename;

	if( filedlg.DoModal() == IDOK)
	{
		filename = filedlg.GetPathName();
                ofstream mystream(filename);

                for (int i=1; i<=pDoc->ndata; i++){
					if (pDoc->x[i]>=Xmin && pDoc->x[i]<=Xmax) {
                    mystream << pDoc->x[i]<< ' ' <<pDoc->y[i] << ' ' << Ycalc[i] << endl;
					}
                }

		mystream.close();
	}

}

void CKawView::OnUpdatebg(CCmdUI* pCmdUI) 
{
	// TODO: この位置に command update UI ハンドラ用のコードを追加してください
	pCmdUI->Enable((F_calc==1) ? TRUE: FALSE);	
	
}