powerchart.cpp

给定一组路面不平度数据

#include <math.h>
#include <stdio.h>

double	realintegrate(double q1,double q0, double x1,double x0,double omiga,double deltl);
double	imagintegrate(double q1,double q0, double x1,double x0,double omiga,double deltl);

void main()
{
	FILE*	fp;
	char	ch[10];
	int		i;
	int		j;
	float	deltl;
	int		totall;
	int		input[1001];
	double	x[1001];
	double	q[1001];
	double	omiga[1000];
	double	sq[1000];
	double	real;
	double	imag;
	fp=fopen("Worse2.txt","r");
	fscanf(fp,"%s%s%f%s%d%s%s%s%",ch,ch,&deltl,ch,&totall,ch,ch,ch);	
	for(i=0;i<=1000;i++)
	{
		fscanf(fp,"%d,",&input[i]);
	}
	fclose(fp);
	for(i=0;i<=1000;i++)
	{
		x[i]=0.25*i;
		q[i]=input[i]/1000.;
	}
	for(j=0;j<=999;j++)
	{
		omiga[j]=2*3.1415926536/deltl/1000*(j+1);
	}
	for(j=0;j<=999;j++)
	{
		real=0.;
		imag=0.;
		for(int i=0;i<=999;i++)
		{
			real+=realintegrate(double(q[i+1]),double(q[i]),x[i+1],x[i],omiga[j],double(deltl));
			imag+=imagintegrate(double(q[i+1]),double(q[i]),x[i+1],x[i],omiga[j],double(deltl));
		}
		sq[j]=(real*real+imag*imag)/(2.*totall*3.1415926536);
	}
	fp=fopen("D:\\project\\matlab\\dat.dat","w");
	for(j=0;j<=999;j++)
	{
		fprintf(fp,"%.12f	%.12f\n",omiga[j],sq[j]);
	}
	fclose(fp);
}
double	realintegrate(double q1,double q0, double x1,double x0,double omiga,double deltl)
{
	double	result;
	double	slope;
	double	factor1;
	double	factor2;
	double	factor3;

	slope=(q1-q0)/deltl;
	factor1=sin(omiga*x1)-sin(omiga*x0);
	factor2=x1*sin(omiga*x1)-x0*sin(omiga*x0);
	factor3=cos(omiga*x1)-cos(omiga*x0);
	result=q0/omiga*factor1-slope/omiga*x0*factor1+slope/omiga*factor2+slope/(omiga*omiga)*factor3;
	return result;
}
double	imagintegrate(double q1,double q0, double x1,double x0,double omiga,double deltl)
{
	double	result;
	double	slope;
	double	factor1;
	double	factor2;
	double	factor3;
	slope=(q1-q0)/deltl;
	factor1=cos(omiga*x0)-cos(omiga*x1);
	factor2=x0*cos(omiga*x0)-x1*cos(omiga*x1);
	factor3=sin(omiga*x1)-sin(omiga*x0);
	result=q0/omiga*factor1-slope/omiga*x0*factor1+slope/omiga*factor2+slope/(omiga*omiga)*factor3;
	return result;
}