#include <iostream>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define PI 3.1415926
#define N 16000 //时间步数
#define oa 0.07
#define ob 0.09
#define intervel 0.05 //时间间隔
#define g 9.8 //重力加速度
#define x 0 //位置
#define h 3 //水深
#define Hs 4.0
#define yy 1.0
#define f0 0.27
#define enlarge 5 //高频相对于谱峰频率的放大倍数
#define mm 1000 //频率间隔
using namespace std;
const int M=N/2;
double S[M],f[M],jiaodu[M],A[M],k[M],nu[N],c[N];
/*****************/
double Jonswap(double f)
{
double JON,o6;
double B=0.0624/(0.23+0.0336*yy-0.185/(1.9+yy));
if(f>f0)o6=ob;
else o6=oa;
JON=B*Hs*Hs*f0*f0*f0*f0/(f*f*f*f*f)*exp(-1.25*(f0*f0*f0*f0)/(f*f*f*f))*pow(yy,exp(-(f-f0)*(f-f0)/(2*o6*o6*f0*f0)));
return(JON);
}
/**********************/
double distributeofw(double a1,double a2,int accuracy)
{
double tt=rand()/(RAND_MAX+1.0);
return (a1+(a2-a1)/accuracy*tt);
}
/*******************/
void cov()
{
for(int i=0;i<N;i++)
{
double sum=0;
for(int j=0;j<N-i;j++)
{
sum=sum+nu[j]*nu[j+i];
}
c[i]=sum/((N-i)*1.0);
}
}
/**********************************/
double disper(double w) //色散方程求不同频率对应的k
{
double answer=w*w/g;
while(abs(answer*tanh(answer*h)/(w*w/g)-1)>0.1e-9)
{
answer=w*w/g/tanh(answer*h);
}
return(answer);
}
/*********************/
void outputvector(FILE *fp,int n,double *vec)
{
for(int i=0;i<n;i++)
{
fprintf(fp,"%.8lf\n",vec[i]);
}
}
/************************/
int main()
{
double Jonswap(double f);
double disper(double w);
double distributeofw(double a1,double a2,int accuracy);
srand( (unsigned)time( NULL ) );
for (int i=0;i<M;i++)
{
double ran=rand()/(RAND_MAX+1.0);
jiaodu[i]=ran*2*PI;//相位角在0~2Π随机分布
}
A[0]=0;f[0]=0;k[0]=0;S[0]=0;
double fH,deltaf;
fH=f0*enlarge*1.0; //右侧频率
deltaf=fH/(1.0*mm); //频率间隔
for (int i=1;i<mm;i++)
{
f[i]=distributeofw((i-1)*deltaf,i*deltaf,100);
k[i]=disper(2*PI*f[i]);
A[i]=sqrt(2*deltaf*Jonswap(f[i]));
}
for(int i=1;i<mm;i++)
{
S[i]=Jonswap(i*deltaf);
}
for(int i=0;i<N;i++)
{
double sum=0;
for(int m=0;m<mm;m++)
{
sum=sum+A[m]*cos(2*PI*f[m]*(i)*intervel-k[m]*x-jiaodu[m]);
}
nu[i]=sum;
}
cov();
/* for(int i=0;i<N;i++)
{
c[i]=0;
for(int j=0;j<N;j++)
{
double fj=1/(N*1.0)/intervel*(j+j+1)/2;
c[i]=c[i]+Jonswap(fj)*cos(2*PI*fj*i*intervel)/(N*1.0)/intervel;
}
}*/
FILE *outputnu;
if((outputnu=fopen("nu.txt","w+"))==NULL)
{
printf("error\n");
exit(0);
}
outputvector(outputnu,N,nu);
fclose(outputnu);
FILE *ouputS;
if((ouputS=fopen("Sf.txt","w+"))==NULL)
{
printf("error\n");
exit(0);
}
outputvector(ouputS,mm,S);
fclose(ouputS);
FILE *ouputc;
if((ouputc=fopen("c.txt","w+"))==NULL)
{
printf("error\n");
exit(0);
}
outputvector(ouputc,N,c);
fclose(ouputc);
return 0;
}
评论0