⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 bysjdoc.cpp

📁 这是关于测井岩心插值的算法
💻 CPP
字号:
// bysjDoc.cpp : implementation of the CBysjDoc class
//

#include "stdafx.h"
#include "bysj.h"

#include "bysjDoc.h"
#include "CntrItem.h"

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

/////////////////////////////////////////////////////////////////////////////
// CBysjDoc

IMPLEMENT_DYNCREATE(CBysjDoc, COleDocument)

BEGIN_MESSAGE_MAP(CBysjDoc, COleDocument)
	//{{AFX_MSG_MAP(CBysjDoc)
	ON_COMMAND(ID_OPEN_DATAFILE, OnOpenDatafile)
	ON_COMMAND(ID_ADD_POINT, OnAddPoint)
	ON_COMMAND(ID_OPEN_ACDATA, OnOpenAcdata)
	ON_COMMAND(ID_FENXI_JIAOZH, OnFenxiJiaozh)
	//}}AFX_MSG_MAP
	// Enable default OLE container implementation
	ON_UPDATE_COMMAND_UI(ID_EDIT_PASTE, COleDocument::OnUpdatePasteMenu)
	ON_UPDATE_COMMAND_UI(ID_EDIT_PASTE_LINK, COleDocument::OnUpdatePasteLinkMenu)
	ON_UPDATE_COMMAND_UI(ID_OLE_EDIT_CONVERT, COleDocument::OnUpdateObjectVerbMenu)
	ON_COMMAND(ID_OLE_EDIT_CONVERT, COleDocument::OnEditConvert)
	ON_UPDATE_COMMAND_UI(ID_OLE_EDIT_LINKS, COleDocument::OnUpdateEditLinksMenu)
	ON_COMMAND(ID_OLE_EDIT_LINKS, COleDocument::OnEditLinks)
	ON_UPDATE_COMMAND_UI_RANGE(ID_OLE_VERB_FIRST, ID_OLE_VERB_LAST, COleDocument::OnUpdateObjectVerbMenu)
	ON_COMMAND(ID_FILE_SEND_MAIL, OnFileSendMail)
	ON_UPDATE_COMMAND_UI(ID_FILE_SEND_MAIL, OnUpdateFileSendMail)
END_MESSAGE_MAP()

BEGIN_DISPATCH_MAP(CBysjDoc, COleDocument)
	//{{AFX_DISPATCH_MAP(CBysjDoc)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//      DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_DISPATCH_MAP
END_DISPATCH_MAP()

// Note: we add support for IID_IBysj to support typesafe binding
//  from VBA.  This IID must match the GUID that is attached to the 
//  dispinterface in the .ODL file.

// {32F70CA4-B831-11D9-A19A-0050BA6DB285}
static const IID IID_IBysj =
{ 0x32f70ca4, 0xb831, 0x11d9, { 0xa1, 0x9a, 0x0, 0x50, 0xba, 0x6d, 0xb2, 0x85 } };

BEGIN_INTERFACE_MAP(CBysjDoc, COleDocument)
	INTERFACE_PART(CBysjDoc, IID_IBysj, Dispatch)
END_INTERFACE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CBysjDoc construction/destruction

CBysjDoc::CBysjDoc()
{
	// Use OLE compound files
	EnableCompoundFile();

	// TODO: add one-time construction code here

	EnableAutomation();

	AfxOleLockApp();
}

CBysjDoc::~CBysjDoc()
{
	AfxOleUnlockApp();
}

BOOL CBysjDoc::OnNewDocument()
{
	if (!COleDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CBysjDoc serialization

void CBysjDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}

	// Calling the base class COleDocument enables serialization
	//  of the container document's COleClientItem objects.
	COleDocument::Serialize(ar);
}

/////////////////////////////////////////////////////////////////////////////
// CBysjDoc diagnostics

#ifdef _DEBUG
void CBysjDoc::AssertValid() const
{
	COleDocument::AssertValid();
}

void CBysjDoc::Dump(CDumpContext& dc) const
{
	COleDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CBysjDoc commands

void CBysjDoc::OnOpenDatafile() 
{
	// TODO: Add your command handler code here
	//打开文件对话框
	CFileDialog dlg(TRUE, ("yxz文件"), NULL,
	OFN_FILEMUSTEXIST | OFN_HIDEREADONLY,
	("yxz文件(*.yxz)|*.yxz|所有文件(*.*)|*.*||"));

	char lpszTitle[]={"读取yxz格式数据文件"};
	dlg.m_ofn.lpstrTitle=lpszTitle;
	if (dlg.DoModal() != IDOK)
	return;    
	m_pathname = dlg.GetPathName();

	if ((fp = fopen(m_pathname,"r"))==NULL)   
	{         
	CString str;
	str="文件"+m_stzFilename+"打开错误";
	AfxMessageBox(str);
	exit(0);
	}
		
	//读数据头
	fseek(fp,43L,0);
	char line[11];
	CString stri  ;
	//把曲线名存入数组
	for(int s=0;s<2;s++)
	{ 
	fgets(line,11,fp);
    stri = line;
	stri.TrimLeft() ; 
    stri.TrimRight();
	linnames[s]=stri;
	}

	//读岩性值数据
	for(int k=0;k<71;k++)
	{
		fscanf(fp,"%f %f",
		&DEPTH1[k],&por[k]);

	}
                                  
	
	rewind(fp);
		fseek(fp,43L,0);
	char lines[10];
	CString strii ;
    	for(int sp=0;sp<2;sp++)
	{ 
	fgets(line,11,fp);
    strii = lines;
	strii.TrimLeft() ; 
    strii.TrimRight();
	linnames[sp]=strii;
	}

	m_Data = TRUE;

	UpdateAllViews(NULL);
}

//画曲线
 void CBysjDoc::draw(CDC*dc)
 {
    //画岩性值曲线	
	float b;
    CString ndata;
	CPen q_pen(PS_SOLID,1,RGB(0,0,255));
	dc->SelectObject(&q_pen);
	for( int l=0;l<(int)70;l++)
	{
		dc->MoveTo(int(90.0+3.0*por[l]),int(90.f+(DEPTH1[l]-1280)*16.0));
		dc->LineTo(int(90.0+3.0*por[l+1]),int(90.f+(DEPTH1[l+1]-1280)*16.0));
	} 
    
	if (INTER == TRUE)
	{
		//插值后的岩性值曲线
		CPen s_pen(PS_SOLID,1,RGB(255,0,255));
		dc->SelectObject(&s_pen);
		for( int t=0;t<318;t++)

		{   
			b=INTER1[t];
			dc->MoveTo(int(90.0+3.f*INTER1[t]),int(90.f+(DD[t]-1280)*16.0));
			dc->LineTo(int(90.0+3.f*INTER1[t+1]),int(90.f+(DD[t+1]-1280)*16.0));
		}					
	}
	//校正后的曲线
	if(m_Data1==TRUE)
	{
	    double  dep=nn*0.125;
		ndata.Format("%f %f %d",dep,max,n1);
			dc->TextOut(630,60,ndata);
			CPen s_pen(PS_SOLID,1,RGB(255,0,255));
		dc->SelectObject(&s_pen);
		for( int t=0;t<309;t++)
		{   
			b=INTER1[t];
			dc->MoveTo(int(500.0+3.f*INTER1[t]),int(90.f+(DD2[t]-1280)*16.0));
			dc->LineTo(int(500.0+3.f*INTER1[t+1]),int(90.f+(DD2[t+1]-1280)*16.0));
		}
			CPen bb_pen(PS_SOLID,1,RGB(0,0,0));
		dc->SelectObject(&bb_pen);
			for(t=0;t<295;t++)
		{   
			b=INTER1[t];
			dc->MoveTo(int(500.0+3.f*yy[t]),int(90.f+(DD2[t]-1280)*16.0));
			dc->LineTo(int(500.0+3.f*yy[t+1]),int(90.f+(DD2[t+1]-1280)*16.0));
		}	


	}
	//测井曲线
	CPen n_pen(PS_SOLID,1,RGB(0,0,255));
	dc->SelectObject(&n_pen);
	for( int ll=0;ll<(int)309;ll++)
	{
		dc->MoveTo(int(10+1.0*AC[ll]),int(90.f+(DEPTH2[ll]-1280)*16.0));
		dc->LineTo(int(10+1.0*AC[ll+1]),int(90.f+(DEPTH2[ll+1]-1280)*16.0));
	} 

 
} 
 
//样条插值程序
void CBysjDoc::OnAddPoint() 

{
	// TODO: Add your command handler code here

	int n,i,j,hh,nn;
	double c[10],d[10],y2[10];
	double g,f,e;
	int klo,khi;
	double h,a,b,aaa,bbb;
	int k=0;
	hh=0;
	
	for(n=1;n<=71;n++)
	{
		for(i=0;i<=3;i++)
		{
			c[i]=(double)DEPTH1[i+n-1];
			d[i]=(double)por[i+n-1];
		}
		spline(c,d,4,0.0,0.0,y2);
		f=y2[0];
		if((int(DEPTH1[n]*8))*0.125==DEPTH1[n])
		{
			nn=int(DEPTH1[n+1]*8)-int(DEPTH1[n]*8)+1;
		}
		else 
		{
            nn=int(DEPTH1[n+1]*8)-int(DEPTH1[n]*8);
		}
	
            for(j=0;j<nn;j++)
			{
				if((int(DEPTH1[n]*8))*0.125==DEPTH1[n])
				{
					g=double(DEPTH1[n]+0.125*j);
				}
				else
				{
					g=double(int((DEPTH1[n]/0.125)+1)*0.125)+0.125*j;
				}
				DD[j+hh]=float(g);
				khi=2;
				klo=1;
				h=c[khi]-c[klo];
				a=(c[khi]-g)/h;
				b=(g-c[klo])/h;
				aaa=a*d[klo]+b*d[khi];
				bbb=(a*a*a-a)*y2[klo]+(b*b*b-b)*y2[khi];
				e=aaa+bbb*h*h/6.0;
				INTER1[j+hh]=(float)e;			
			}	
		hh=hh+nn;
	}
    FILE *fp5;
	if((fp5=fopen("E:\\gang\\bysj\\ch25cz.txt","w"))==NULL)
	exit(0);
    for(j=0;j<309;j++)
	{
		fprintf(fp5,"%7.3f  %7.3f\n",DD[j],INTER1[j]);
	}
	fclose(fp5);
	INTER = TRUE;
	UpdateAllViews(NULL);
}

/////////////////////////////////////////////////////////////////////
void CBysjDoc::spline(double x[],double y[],int n,double yp1,double ypn,double y2[])
{
	double u[1001],aaa,sig,p,bbb,ccc,qn,un;
    int i,k;
	if(yp1>9.9e+29)
	{
		y2[0]=0;
		u[0]=0;
	}
	else
	{
		y2[0]=-0.5;
		aaa=(y[1]-y[0])/(x[1]-x[0]);
		u[0]=(3.0/(x[1]-x[0]))*(aaa-yp1);
	}
	for(i=1;i<=n-2;i++)
	{
		sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
		p=sig*y2[i-1]+2.0;
		y2[i]=(sig-1.0)/p;
		aaa=(y[i+1]-y[i])/(x[i+1]-x[i]);
		bbb=(y[i]-y[i-1])/(x[i]-x[i-1]);
		ccc=x[i+1]-x[i-1];
		u[i]=(6.0*(aaa-bbb)/ccc-sig*u[i-1])/p;
	}
	if(ypn>9.9e+29)
	{
		qn=0.0;
		un=0.0;
	}
	else
	{
		qn=0.5;
		aaa=ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
		un=(3.0/(x[n-1]-x[n-2]))*aaa;
	}
	y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
	for(k=n-2;k>=0;k--)
	{
		y2[k]=y2[k]*y2[k+1]+u[k];
	}

}
//读入测井数据
void CBysjDoc::OnOpenAcdata() 
{
	// TODO: Add your command handler code here
		//打开文件对话框
	CFileDialog dlg(TRUE, ("txt文件"), NULL,
	OFN_FILEMUSTEXIST | OFN_HIDEREADONLY,
	("txt文件(*.txt)|*.txt|所有文件(*.*)|*.*||"));

	char lpszTitle[]={"读取txt格式数据文件"};
	dlg.m_ofn.lpstrTitle=lpszTitle;
	if (dlg.DoModal() != IDOK)
	return;    
	m_pathname = dlg.GetPathName();

	if ((fp1= fopen(m_pathname,"r"))==NULL)   
	{         
	CString str;
	str="文件"+m_stzFilename+"打开错误";
	AfxMessageBox(str);
	exit(0);
	}		
	//读数据头
	fseek(fp1,43L,0);
	char line[10];
	CString stri;
	//把曲线名存入数组
	for(int s=0;s<3;s++)
	{ 
	fgets(line,10,fp1);
    stri = line;
	stri.TrimLeft() ;       
    stri.TrimRight();
	linnames1[s]=stri;
	}
	//读取测井值
	for(int k=0;k<310;k++)
	{
        fscanf(fp1,"%f %f %f",
		&DEPTH2[k],&SP[k],&AC[k]);
	}                          
	rewind(fp1);
	m_Data = TRUE;
	UpdateAllViews(NULL);
}
//相关分析深度校正
void CBysjDoc::OnFenxiJiaozh() 
{
	// TODO: Add your command handler code here
	double MY[120],VY[120],SY[120],CXY[120],C[120];
	double MX,VX,SX;
	double a,b,c,d,e;
	int    kk,ll;
	a=0.0;
	kk=40;
	ll=200;
	int mm=ll-kk;
	for(int m=kk;m<ll;m++)
	{
		a=a+AC[m];
	}
	MX=a/float(mm);
	c=0.0;
	for(m=kk;m<ll;m++)
	{
		c=c+(AC[m]-MX)*(AC[m]-MX);
	}
	VX=c/float(mm);
	SX=sqrt(VX);
	for(int n=0;n<120;n++)
	{
		b=0.0;
		for(int i=0;i<mm;i++)
		{
			b=b+INTER1[n+i];
		}
		MY[n]=b/float(mm);
		d=0.0;
	    for(i=0;i<mm;i++)
		{
			d=d+(INTER1[n+i]-MY[n])*(INTER1[n+i]-MY[n]);
		}                                                  
		VY[n]=d/float(mm);
		SY[n]=sqrt(VY[n]);
	    e=0.0;
		for(i=0;i<mm;i++)
		{
		   e=e+(AC[i+kk]-MX)*(INTER1[n+i]-MY[n]);
		}
		CXY[n]=e/float(mm);
		C[n]=CXY[n]/(SX*SY[n]);
		if(C[n]>0.0)
		{
			C[n]=C[n];
		}
		else
		{
			C[n]=(-1.0)*C[n];
		}
	}
	 max=C[0];
	for(n=0;n<120;n++)
	{
		if(C[n]>max)
		{
			max=C[n];
			nn=n;
		}
	}
	int qq=int((DD[0]-DEPTH2[0])/0.125);
	nn=kk-nn-qq;
	for(int ii=0;ii<310;ii++)
	{
		DD2[ii]=DD[ii]+float(0.125*nn);
	}	
	m_Data1 = TRUE;
	UpdateAllViews(NULL);
	float aa=DD2[0];
	for(ii=0;ii<310;ii++)
	{
		if(aa==DEPTH2[ii])
		{
			n1=ii;
		}
	}
/*	double x=0.0;
	for(ii=n1;ii<309;ii++)
	{
		x=x+AC[ii];
	}
	x=x/double(309-n1);
	double y=0.0;
	int ab=309-n1;
	for(ii=60;ii<ab;ii++)
	{
		y=y+INTER1[ii];
	}
	y=y/double(309-n1);
    double b2=0.0;
    double b3=0.0;
	for(ii=0;ii<(309-n1);ii++)
	{
		b2=b2+(AC[ii+n1]-x)*(INTER1[ii]-y);
		b3=b3+(AC[ii+n1]-x)*(AC[ii+n1]-x);
	}
	b1=b2/b3;
	b0=y-b1*x;*/

	double wucha[300];
	for(ii=0;ii<(309-n1);ii++)
	{
		yy[ii]=float(0.186*AC[n1+ii]-37.765);
		wucha[ii]=yy[ii]-INTER1[ii];
	}
	for(ii=0;ii<(309-n1);ii++)
	{
		if(wucha[ii]<0)
		{
			wucha[ii]=(-1)*wucha[ii];
		}
	}
	FILE *fp3;
	if((fp3=fopen("E:\\GILCEP\\Bin\\yinter.txt","w"))==NULL)
	exit(0);
	for(ii=60;ii<210;ii++)
	{
		if(wucha[ii]<5.0)
		{
			fprintf(fp3,"%7.3f  %7.3f  %7.3f\n",DEPTH2[n1+ii],yy[ii],INTER1[ii]);
		}
	}
    fclose(fp3);
	//输出建模文件
    FILE *fp2;
	if((fp2=fopen("E:\\GILCEP\\Bin\\acp.txt","w"))==NULL)
	exit(0);
	for(ii=60;ii<210;ii++)
	{
		if(INTER1[ii]>0&&(AC[ii+n1]<345.0||INTER1[ii]>16.0))
		{
			fprintf(fp2,"%7.3f  %7.3f  %7.3f\n",DEPTH2[ii+n1],INTER1[ii],AC[ii+n1]);
		}
	}
	fclose(fp2);

	//误差分析
	float ff=0.0;
	int ss,dd;
	ff=0.0;
	ss=0;
	//DEPTN1(1300,1302)
    for(ii=10;ii<=24;ii++)
	{   
			ff=ff+por[ii];
			ss=ss+1;	
	}
	ff=ff/float(ss);

	float ff1=0.0;
	ss=0;
	//DEPTN1(1307,1311)
    for(ii=27;ii<=43;ii++)
	{   
			ff1=ff1+por[ii];
			ss=ss+1;	
	}
	ff1=ff1/float(ss);

	//DEPTN2(1300.5,1302.5)
	float fff=0.0;
    dd=0;
	for(ii=104;ii<120;ii++)
	{		
		fff=fff+yy[ii];		
	     dd=dd+1;
	}
	fff=fff/float(dd);

	//DEPTN2(1307.5,1311.5)
	float fff1=0.0;
    dd=0;
	for(ii=160;ii<=192;ii++)
	{		
		fff1=fff1+yy[ii];		
	      dd=dd+1;
	}
	fff1=fff1/float(dd);
	
	float cha=ff-fff;
    float cha1=ff1-fff1;
    FILE *fp4;
	if((fp4=fopen("E:\\GILCEP\\Bin\\wucha.txt","w"))==NULL)
		exit(0);
	fprintf(fp4,"1300-1302 %7.3f  %7.3f  %7.3f\n1307-1311 %7.3f  %7.3f  %7.3f\n",ff,fff,cha,ff1,fff1,cha1);
    fclose(fp4);
		//输出建模文件
    //1293-1301
    FILE *fp6;
	if((fp6=fopen("E:\\GILCEP\\Bin\\acp1.txt","w"))==NULL)
	exit(0);
	for(ii=44;ii<108;ii++)
	{	
		fprintf(fp6,"%7.3f  %7.3f  %7.3f\n",DEPTH2[ii+n1],INTER1[ii],yy[ii]);
	}
	fclose(fp6);
    //1303-1307
    FILE *fp7;
	if((fp7=fopen("E:\\GILCEP\\Bin\\acp2.txt","w"))==NULL)
	exit(0);
	for(ii=124;ii<156;ii++)    
	{
	
		fprintf(fp7,"%7.3f  %7.3f  %7.3f\n",DEPTH2[ii+n1],INTER1[ii],yy[ii]);
	
	}
	fclose(fp7);
    //1293-1307
    FILE *fp8;
	if((fp8=fopen("E:\\GILCEP\\Bin\\acp3.txt","w"))==NULL)
	exit(0);
	for(ii=44;ii<156;ii++)    
	{
	
		fprintf(fp8,"%7.3f  %7.3f  %7.3f\n",DEPTH2[ii+n1],INTER1[ii],yy[ii]);
	
	}
	fclose(fp8);
}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -