danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Restored c_functions.c

Browse Source
master
Per A Brodtkorb 9 years ago
parent cd7b75b184
commit 1863b2233f
1 changed files with 313 additions and 311 deletions
Show Stats Download Patch File Download Diff File

624
wafo/source/c_library/c_functions.c
Unescape Escape View File

@ -1,11 +1,11 @@
#include "math.h" #include "math.h"
/* /*
* Install gfortran and run the following to build the module on windows: * Install gfortran and run the following to build the module on windows:
* f2py c_library.pyf c_functions.c -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 * f2py c_library.pyf c_functions.c -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
*/ */
/* /*
* findrfc.c - * findrfc.c -
* *
* Returns indices to RFC turningpoints of a vector * Returns indices to RFC turningpoints of a vector
* of turningpoints * of turningpoints
@ -15,11 +15,11 @@
void findrfc(double *y1,double hmin, int *ind, int n,int *info) { void findrfc(double *y1,double hmin, int *ind, int n,int *info) {
double xminus,xplus,Tpl,Tmi,*y,Tstart; double xminus,xplus,Tpl,Tmi,*y,Tstart;
int i,j,ix=0,NC; int i,j,ix=0,NC,iy;
info[0] = 0; info[0] = 0;
if (*(y1+0)> *(y1+1)){ if (*(y1+0)> *(y1+1)){
/* if first is a max , ignore the first max*/ /* if first is a max , ignore the first max*/
y=&(*(y1+1)); y=&(*(y1+1));
NC=floor((n-1)/2); NC=floor((n-1)/2);
Tstart=1; Tstart=1;
} }
@ -28,140 +28,142 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) {
NC=floor(n/2); NC=floor(n/2);
Tstart=0; Tstart=0;
} }
if (NC<1){ if (NC<1){
return; /* No RFC cycles*/ return; /* No RFC cycles*/
} }
if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){ if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){
info[0] = -1; info[0] = -1;
return; /*This is not a sequence of turningpoints, exit */ return; /*This is not a sequence of turningpoints, exit */
} }
if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){ if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){
info[0]=-1; info[0]=-1;
return; /*This is not a sequence of turningpoints, exit */ return; /*This is not a sequence of turningpoints, exit */
} }
for (i=0; i<NC; i++) { for (i=0; i<NC; i++) {
Tmi=Tstart+2*i; Tmi=Tstart+2*i;
Tpl=Tstart+2*i+2; Tpl=Tstart+2*i+2;
xminus=*(y+2*i); xminus=*(y+2*i);
xplus=*(y+2*i+2); xplus=*(y+2*i+2);
if(i!=0){ if(i!=0){
j=i-1; j=i-1;
while((j>=0) && (*(y+2*j+1)<=*(y+2*i+1))){ while((j>=0) && (*(y+2*j+1)<=*(y+2*i+1))){
if( (*(y+2*j)<xminus) ){ if( (*(y+2*j)<xminus) ){
xminus=*(y+2*j); xminus=*(y+2*j);
Tmi=Tstart+2*j; Tmi=Tstart+2*j;
} /*if */ } /*if */
j--; j--;
} /*while j*/ } /*while j*/
} /*if i */ } /*if i */
if ( xminus >= xplus){ if ( xminus >= xplus){
if ( (*(y+2*i+1)-xminus) >= hmin){ if ( (*(y+2*i+1)-xminus) >= hmin){
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
goto L180; goto L180;
} }
j=i+1; j=i+1;
while((j<NC) ) { while((j<NC) ) {
if (*(y+2*j+1) >= *(y+2*i+1)) goto L170; if (*(y+2*j+1) >= *(y+2*i+1)) goto L170;
if( (*(y+2*j+2) <= xplus) ){ if( (*(y+2*j+2) <= xplus) ){
xplus=*(y+2*j+2); xplus=*(y+2*j+2);
Tpl=(Tstart+2*j+2); Tpl=(Tstart+2*j+2);
}/*if*/ }/*if*/
j++; j++;
} /*while*/ } /*while*/
if ( (*(y+2*i+1)-xminus) >= hmin) { if ( (*(y+2*i+1)-xminus) >= hmin) {
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
goto L180; goto L180;
L170: L170:
if (xplus <= xminus ) { if (xplus <= xminus ) {
if ( (*(y+2*i+1)-xminus) >= hmin){ if ( (*(y+2*i+1)-xminus) >= hmin){
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
/*goto L180;*/ /*goto L180;*/
} }
else{ else{
if ( (*(y+2*i+1)-xplus) >= hmin) { if ( (*(y+2*i+1)-xplus) >= hmin) {
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
*(ind+ix)=Tpl; *(ind+ix)=Tpl;
ix++; ix++;
} /*if*/ } /*if*/
} /*elseif*/ } /*elseif*/
L180: L180:
iy=i;
} /* for i */ } /* for i */
info[0] = ix; info[0] = ix;
return ; return ;
} }
/* /*
* findcross.c - * findcross.c -
* *
* Returns indices to level v crossings of argument vector * Returns indices to level v crossings of argument vector
* *
* 1998 by Per Andreas Brodtkorb. last modified 23.06-98 * 1998 by Per Andreas Brodtkorb. last modified 23.06-98
*/ */
void findcross(double *y, double v, int *ind, int n, int *info) void findcross(double *y, double v, int *ind, int n, int *info)
{ int i,start, ix=0,dcross=0; { int i,start, ix=0,dcross=0;
start=0; start=0;
if ( y[0]< v){ if ( y[0]< v){
dcross=-1; /* first is a up-crossing*/ dcross=-1; /* first is a up-crossing*/
} }
else if ( y[0]> v){ else if ( y[0]> v){
dcross=1; /* first is a down-crossing*/ dcross=1; /* first is a down-crossing*/
} }
else if ( y[0]== v){ else if ( y[0]== v){
/* Find out what type of crossing we have next time.. */ /* Find out what type of crossing we have next time.. */
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
start=i; start=i;
if ( y[i]< v){ if ( y[i]< v){
ind[ix] = i-1; /* first crossing is a down crossing*/ ind[ix] = i-1; /* first crossing is a down crossing*/
ix++; ix++;
dcross=-1; /* The next crossing is a up-crossing*/ dcross=-1; /* The next crossing is a up-crossing*/
goto L120; goto L120;
} }
else if ( y[i]> v){ else if ( y[i]> v){
ind[ix] = i-1; /* first crossing is a up-crossing*/ ind[ix] = i-1; /* first crossing is a up-crossing*/
ix++; ix++;
dcross=1; /*The next crossing is a down-crossing*/ dcross=1; /*The next crossing is a down-crossing*/
goto L120; goto L120;
} }
} }
} }
L120: L120:
for (i=start; i<n-1; i++) { for (i=start; i<n-1; i++) {
if (( (dcross==-1) && (y[i]<=v) && (y[i+1] > v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { if (( (dcross==-1) && (y[i]<=v) && (y[i+1] > v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) {
ind[ix] = i; ind[ix] = i;
ix++; ix++;
dcross=-dcross; dcross=-dcross;
} }
} }
info[0] = ix; info[0] = ix;
return; return;
} }
@ -171,19 +173,19 @@ void findcross(double *y, double v, int *ind, int n, int *info)
* CALL: disufq(rvec,ivec,rA,iA, w,kw,h,g,nmin,nmax,m,n) * CALL: disufq(rvec,ivec,rA,iA, w,kw,h,g,nmin,nmax,m,n)
* *
* rvec, ivec = real and imaginary parts of the resultant (size m X n). * rvec, ivec = real and imaginary parts of the resultant (size m X n).
* rA, iA = real and imaginary parts of the amplitudes (size m X n). * rA, iA = real and imaginary parts of the amplitudes (size m X n).
* w = vector with angular frequencies (w>=0) * w = vector with angular frequencies (w>=0)
* kw = vector with wavenumbers (kw>=0) * kw = vector with wavenumbers (kw>=0)
* h = water depth (h >=0) * h = water depth (h >=0)
* g = constant acceleration of gravity * g = constant acceleration of gravity
* nmin = minimum index where rA(:,nmin) and iA(:,nmin) is * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is
* greater than zero. * greater than zero.
* nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is
* greater than zero. * greater than zero.
* m = size(rA,1),size(iA,1) * m = size(rA,1),size(iA,1)
* n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2) * n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2)
* *
* DISUFQ returns the summation of difference frequency and sum * DISUFQ returns the summation of difference frequency and sum
* frequency effects in the vector vec = rvec +sqrt(-1)*ivec. * frequency effects in the vector vec = rvec +sqrt(-1)*ivec.
* The 2'nd order contribution to the Stokes wave is then calculated by * The 2'nd order contribution to the Stokes wave is then calculated by
* a simple 1D Fourier transform, real(FFT(vec)). * a simple 1D Fourier transform, real(FFT(vec)).
@ -193,15 +195,15 @@ void findcross(double *y, double v, int *ind, int n, int *info)
* *
* by Per Andreas Brodtkorb 15.08.2001 * by Per Andreas Brodtkorb 15.08.2001
* revised pab 14.03.2002, 01.05.2002 22.07.2002, oct 2008 * revised pab 14.03.2002, 01.05.2002 22.07.2002, oct 2008
*/ */
void disufq(double *rvec, double *ivec, void disufq(double *rvec, double *ivec,
double *rA, double *iA, double *rA, double *iA,
double *w, double *kw, double *w, double *kw,
double h, double g, double h, double g,
int nmin, int nmax, int nmin, int nmax,
int m, int n) int m, int n)
{ {
double Epij, Edij; double Epij, Edij;
double tmp1, tmp2, tmp3, tmp4, kfact; double tmp1, tmp2, tmp3, tmp4, kfact;
double w1, w2, kw1, kw2, Cg; double w1, w2, kw1, kw2, Cg;
@ -215,7 +217,7 @@ void disufq(double *rvec, double *ivec,
} }
// kfact is set to 2 in order to exploit the symmetry. // kfact is set to 2 in order to exploit the symmetry.
// If you set kfact to 1, you must uncomment all statements // If you set kfact to 1, you must uncomment all statements
// including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2]. // including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
kfact = 2.0; kfact = 2.0;
@ -225,70 +227,70 @@ void disufq(double *rvec, double *ivec,
iz1 = 2*ixi; iz1 = 2*ixi;
//iz2 = n*m-ixi; //iz2 = n*m-ixi;
kw1 = kw[ix]; kw1 = kw[ix];
Epij = kw1; Epij = kw1;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
tmp1 = kfact*(rrA-iiA)*Epij; tmp1 = kfact*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*riA*Epij; tmp2 = kfact*2.0*riA*Epij;
rvec[iz1] += tmp1; rvec[iz1] += tmp1;
ivec[iz1] += tmp2; ivec[iz1] += tmp2;
//rvec[iz2] += tmp1; //rvec[iz2] += tmp1;
//ivec[iz2] -= tmp2; //ivec[iz2] -= tmp2;
//iz2++; //iz2++;
/// Difference frequency effects are zero along the diagonal /// Difference frequency effects are zero along the diagonal
/// and are thus not contributing to the mean. /// and are thus not contributing to the mean.
} }
for (jy = ix+1;jy<nmax;jy++){ for (jy = ix+1;jy<nmax;jy++){
kw2 = kw[jy]; kw2 = kw[jy];
Epij = 0.5*(kw2 + kw1); Epij = 0.5*(kw2 + kw1);
Edij = -0.5*(kw2 - kw1); Edij = -0.5*(kw2 - kw1);
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
jyi = jy*m; jyi = jy*m;
iz1 = ixi+jyi; iz1 = ixi+jyi;
iv1 = jyi-ixi; iv1 = jyi-ixi;
//iz2 = (n*m-iz1); //iz2 = (n*m-iz1);
//iv2 = (n*m-iv1); //iv2 = (n*m-iv1);
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1; rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2; ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
//rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1; //rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
//ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2; //ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2;
// iz2++; // iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2; ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
//rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1; //rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
//ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2; //ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
//iv2++; //iv2++;
} }
} }
} }
} }
else{ /* Finite water depth */ else{ /* Finite water depth */
for (ix = nmin-1;ix<nmax;ix++) { for (ix = nmin-1;ix<nmax;ix++) {
kw1 = kw[ix]; kw1 = kw[ix];
w1 = w[ix]; w1 = w[ix];
tmp1 = tanh(kw1*h); tmp1 = tanh(kw1*h);
/// Cg, wave group velocity /// Cg, wave group velocity
@ -296,7 +298,7 @@ void disufq(double *rvec, double *ivec,
tmp1 = 0.5*g*(kw1/w1)*(kw1/w1); tmp1 = 0.5*g*(kw1/w1)*(kw1/w1);
tmp2 = 0.5*w1*w1/g; tmp2 = 0.5*w1*w1/g;
tmp3 = g*kw1/(w1*Cg); tmp3 = g*kw1/(w1*Cg);
if (kw1*h<300.0){ if (kw1*h<300.0){
tmp4 = kw1/sinh(2.0*kw1*h); tmp4 = kw1/sinh(2.0*kw1*h);
} }
@ -305,35 +307,35 @@ void disufq(double *rvec, double *ivec,
} }
// Difference frequency effects finite water depth // Difference frequency effects finite water depth
Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK
// Sum frequency effects finite water depth // Sum frequency effects finite water depth
Epij = (3.0*(tmp1-tmp2)/(1.0-tmp1/kw1*tanh(2.0*kw1*h))+3.0*tmp2-tmp1); /// OK Epij = (3.0*(tmp1-tmp2)/(1.0-tmp1/kw1*tanh(2.0*kw1*h))+3.0*tmp2-tmp1); /// OK
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
iz1 = 2*ixi; iz1 = 2*ixi;
//iz2 = n*m-ixi; //iz2 = n*m-ixi;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
rvec[iz1] += kfact*(rrA-iiA)*Epij; rvec[iz1] += kfact*(rrA-iiA)*Epij;
ivec[iz1] += kfact*2.0*riA*Epij; ivec[iz1] += kfact*2.0*riA*Epij;
//rvec[iz2] += kfact*(rrA-iiA)*Epij; //rvec[iz2] += kfact*(rrA-iiA)*Epij;
//ivec[iz2] -= kfact*2.0*riA*Epij; //ivec[iz2] -= kfact*2.0*riA*Epij;
//iz2++; //iz2++;
/// Difference frequency effects along the diagonal /// Difference frequency effects along the diagonal
/// are only contributing to the mean /// are only contributing to the mean
rvec[i] += 2.0*(rrA+iiA)*Edij; rvec[i] += 2.0*(rrA+iiA)*Edij;
} }
for (jy = ix+1;jy<nmax;jy++) { for (jy = ix+1;jy<nmax;jy++) {
// w1 = w[ix]; // w1 = w[ix];
// kw1 = kw[ix]; // kw1 = kw[ix];
w2 = w[jy]; w2 = w[jy];
kw2 = kw[jy]; kw2 = kw[jy];
tmp1 = g*(kw1/w1)*(kw2/w2); tmp1 = g*(kw1/w1)*(kw2/w2);
@ -341,7 +343,7 @@ void disufq(double *rvec, double *ivec,
tmp3 = 0.5*g*(w1*kw2*kw2+w2*kw1*kw1)/(w1*w2*(w1+w2)); tmp3 = 0.5*g*(w1*kw2*kw2+w2*kw1*kw1)/(w1*w2*(w1+w2));
tmp4 = (1-g*(kw1+kw2)/(w1+w2)/(w1+w2)*tanh((kw1+kw2)*h)); tmp4 = (1-g*(kw1+kw2)/(w1+w2)/(w1+w2)*tanh((kw1+kw2)*h));
Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */ Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */
tmp2 = 0.5/g*(w1*w1+w2*w2-w1*w2); /*OK*/ tmp2 = 0.5/g*(w1*w1+w2*w2-w1*w2); /*OK*/
tmp3 = -0.5*g*(w1*kw2*kw2-w2*kw1*kw1)/(w1*w2*(w1-w2)); tmp3 = -0.5*g*(w1*kw2*kw2-w2*kw1*kw1)/(w1*w2*(w1-w2));
tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h)); tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h));
@ -354,30 +356,30 @@ void disufq(double *rvec, double *ivec,
iv1 = jyi-ixi; iv1 = jyi-ixi;
// iz2 = (n*m-iz1); // iz2 = (n*m-iz1);
// iv2 = n*m-iv1; // iv2 = n*m-iv1;
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rvec[iz1] += tmp1;///rvec[i][jy+ix] += tmp1; rvec[iz1] += tmp1;///rvec[i][jy+ix] += tmp1;
ivec[iz1] += tmp2;///ivec[i][jy+ix] += tmp2; ivec[iz1] += tmp2;///ivec[i][jy+ix] += tmp2;
//rvec[iz2] += tmp1;///rvec[i][n*m-(jy+ix)] += tmp1; //rvec[iz2] += tmp1;///rvec[i][n*m-(jy+ix)] += tmp1;
//ivec[iz2] -= tmp2;///ivec[i][n*m-(jy+ix)] -= tmp2; //ivec[iz2] -= tmp2;///ivec[i][n*m-(jy+ix)] -= tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
ivec[iv1] += tmp2;///ivec[i][jy-ix] -= tmp2; ivec[iv1] += tmp2;///ivec[i][jy-ix] -= tmp2;
//rvec[iv2] += tmp1; //rvec[iv2] += tmp1;
//ivec[iv2] -= tmp2; //ivec[iv2] -= tmp2;
//iv2++; //iv2++;
} }
} }
} }
@ -393,19 +395,19 @@ void disufq(double *rvec, double *ivec,
* effects (size m X n). * effects (size m X n).
* rdvec, idvec = real and imaginary parts of the difference frequency * rdvec, idvec = real and imaginary parts of the difference frequency
* effects (size m X n). * effects (size m X n).
* rA, iA = real and imaginary parts of the amplitudes (size m X n). * rA, iA = real and imaginary parts of the amplitudes (size m X n).
* w = vector with angular frequencies (w>=0) * w = vector with angular frequencies (w>=0)
* kw = vector with wavenumbers (kw>=0) * kw = vector with wavenumbers (kw>=0)
* h = water depth (h >=0) * h = water depth (h >=0)
* g = constant acceleration of gravity * g = constant acceleration of gravity
* nmin = minimum index where rA(:,nmin) and iA(:,nmin) is * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is
* greater than zero. * greater than zero.
* nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is
* greater than zero. * greater than zero.
* m = size(rA,1),size(iA,1) * m = size(rA,1),size(iA,1)
* n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2) * n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2)
* *
* DISUFQ2 returns the summation of sum and difference frequency * DISUFQ2 returns the summation of sum and difference frequency
* frequency effects in the vectors svec = rsvec +sqrt(-1)*isvec and * frequency effects in the vectors svec = rsvec +sqrt(-1)*isvec and
* dvec = rdvec +sqrt(-1)*idvec. * dvec = rdvec +sqrt(-1)*idvec.
* The 2'nd order contribution to the Stokes wave is then calculated by * The 2'nd order contribution to the Stokes wave is then calculated by
@ -414,17 +416,17 @@ void disufq(double *rvec, double *ivec,
* *
* This is a MEX-file for MATLAB. * This is a MEX-file for MATLAB.
* by Per Andreas Brodtkorb 15.08.2001 * by Per Andreas Brodtkorb 15.08.2001
* revised pab 14.03.2002, 01.05.2002 * revised pab 14.03.2002, 01.05.2002
*/ */
void disufq2(double *rsvec, double *isvec, void disufq2(double *rsvec, double *isvec,
double *rdvec, double *idvec, double *rdvec, double *idvec,
double *rA, double *iA, double *rA, double *iA,
double *w, double *kw, double *w, double *kw,
double h, double g, double h, double g,
int nmin, int nmax, int nmin, int nmax,
int m, int n) int m, int n)
{ {
double Epij, Edij; double Epij, Edij;
double tmp1, tmp2, tmp3, tmp4, kfact; double tmp1, tmp2, tmp3, tmp4, kfact;
double w1, w2, kw1, kw2, Cg; double w1, w2, kw1, kw2, Cg;
@ -441,7 +443,7 @@ void disufq2(double *rsvec, double *isvec,
} }
// kfact is set to 2 in order to exploit the symmetry. // kfact is set to 2 in order to exploit the symmetry.
// If you set kfact to 1, you must uncomment all statements // If you set kfact to 1, you must uncomment all statements
// including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2]. // including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
kfact = 2.0; kfact = 2.0;
@ -451,70 +453,70 @@ void disufq2(double *rsvec, double *isvec,
iz1 = 2*ixi; iz1 = 2*ixi;
//iz2 = n*m-ixi; //iz2 = n*m-ixi;
kw1 = kw[ix]; kw1 = kw[ix];
Epij = kw1; Epij = kw1;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
tmp1 = kfact*(rrA-iiA)*Epij; tmp1 = kfact*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*riA*Epij; tmp2 = kfact*2.0*riA*Epij;
rsvec[iz1] += tmp1; rsvec[iz1] += tmp1;
isvec[iz1] += tmp2; isvec[iz1] += tmp2;
//rsvec[iz2] += tmp1; //rsvec[iz2] += tmp1;
//isvec[iz2] -= tmp2; //isvec[iz2] -= tmp2;
//iz2++; //iz2++;
/// Difference frequency effects are zero along the diagonal /// Difference frequency effects are zero along the diagonal
/// and are thus not contributing to the mean. /// and are thus not contributing to the mean.
} }
for (jy = ix+1;jy<nmax;jy++){ for (jy = ix+1;jy<nmax;jy++){
kw2 = kw[jy]; kw2 = kw[jy];
Epij = 0.5*(kw2 + kw1); Epij = 0.5*(kw2 + kw1);
Edij = -0.5*(kw2 - kw1); Edij = -0.5*(kw2 - kw1);
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
jyi = jy*m; jyi = jy*m;
iz1 = ixi+jyi; iz1 = ixi+jyi;
iv1 = jyi-ixi; iv1 = jyi-ixi;
//iz2 = (n*m-iz1); //iz2 = (n*m-iz1);
//iv2 = (n*m-iv1); //iv2 = (n*m-iv1);
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rsvec[iz1] += tmp1; ///rvec[i][ix+jy] += tmp1; rsvec[iz1] += tmp1; ///rvec[i][ix+jy] += tmp1;
isvec[iz1] += tmp2; ///ivec[i][ix+jy] += tmp2; isvec[iz1] += tmp2; ///ivec[i][ix+jy] += tmp2;
//rsvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1; //rsvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
//isvec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] += tmp2; //isvec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] += tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rdvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rdvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
idvec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2; idvec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
//rdvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1; //rdvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
//idvec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2; //idvec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
// iv2++; // iv2++;
} }
} }
} }
} }
else{ /* Finite water depth */ else{ /* Finite water depth */
for (ix = nmin-1;ix<nmax;ix++) { for (ix = nmin-1;ix<nmax;ix++) {
kw1 = kw[ix]; kw1 = kw[ix];
w1 = w[ix]; w1 = w[ix];
tmp1 = tanh(kw1*h); tmp1 = tanh(kw1*h);
/// Cg, wave group velocity /// Cg, wave group velocity
@ -522,7 +524,7 @@ void disufq2(double *rsvec, double *isvec,
tmp1 = 0.5*g*(kw1/w1)*(kw1/w1); tmp1 = 0.5*g*(kw1/w1)*(kw1/w1);
tmp2 = 0.5*w1*w1/g; tmp2 = 0.5*w1*w1/g;
tmp3 = g*kw1/(w1*Cg); tmp3 = g*kw1/(w1*Cg);
if (kw1*h<300.0){ if (kw1*h<300.0){
tmp4 = kw1/sinh(2.0*kw1*h); tmp4 = kw1/sinh(2.0*kw1*h);
} }
@ -531,27 +533,27 @@ void disufq2(double *rsvec, double *isvec,
} }
// Difference frequency effects finite water depth // Difference frequency effects finite water depth
Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK
// Sum frequency effects finite water depth // Sum frequency effects finite water depth
Epij = (3.0*(tmp1-tmp2)/(1.0-tmp1/kw1*tanh(2.0*kw1*h))+3.0*tmp2-tmp1); /// OK Epij = (3.0*(tmp1-tmp2)/(1.0-tmp1/kw1*tanh(2.0*kw1*h))+3.0*tmp2-tmp1); /// OK
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
iz1 = 2*ixi; iz1 = 2*ixi;
//iz2 = n*m-ixi; //iz2 = n*m-ixi;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
rsvec[iz1] += kfact*(rrA-iiA)*Epij; rsvec[iz1] += kfact*(rrA-iiA)*Epij;
isvec[iz1] += kfact*2.0*riA*Epij; isvec[iz1] += kfact*2.0*riA*Epij;
//rsvec[iz2] += kfact*(rrA-iiA)*Epij; //rsvec[iz2] += kfact*(rrA-iiA)*Epij;
//isvec[iz2] -= kfact*2.0*riA*Epij; //isvec[iz2] -= kfact*2.0*riA*Epij;
/// Difference frequency effects along the diagonal /// Difference frequency effects along the diagonal
/// are only contributing to the mean /// are only contributing to the mean
//printf(" %f \n",2.0*(rrA+iiA)*Edij); //printf(" %f \n",2.0*(rrA+iiA)*Edij);
@ -559,7 +561,7 @@ void disufq2(double *rsvec, double *isvec,
} }
for (jy = ix+1;jy<nmax;jy++) { for (jy = ix+1;jy<nmax;jy++) {
// w1 = w[ix]; // w1 = w[ix];
// kw1 = kw[ix]; // kw1 = kw[ix];
w2 = w[jy]; w2 = w[jy];
kw2 = kw[jy]; kw2 = kw[jy];
tmp1 = g*(kw1/w1)*(kw2/w2); tmp1 = g*(kw1/w1)*(kw2/w2);
@ -567,7 +569,7 @@ void disufq2(double *rsvec, double *isvec,
tmp3 = 0.5*g*(w1*kw2*kw2+w2*kw1*kw1)/(w1*w2*(w1+w2)); tmp3 = 0.5*g*(w1*kw2*kw2+w2*kw1*kw1)/(w1*w2*(w1+w2));
tmp4 = (1-g*(kw1+kw2)/(w1+w2)/(w1+w2)*tanh((kw1+kw2)*h)); tmp4 = (1-g*(kw1+kw2)/(w1+w2)/(w1+w2)*tanh((kw1+kw2)*h));
Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */ Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */
tmp2 = 0.5/g*(w1*w1+w2*w2-w1*w2); /*OK*/ tmp2 = 0.5/g*(w1*w1+w2*w2-w1*w2); /*OK*/
tmp3 = -0.5*g*(w1*kw2*kw2-w2*kw1*kw1)/(w1*w2*(w1-w2)); tmp3 = -0.5*g*(w1*kw2*kw2-w2*kw1*kw1)/(w1*w2*(w1-w2));
tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h)); tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h));
@ -580,30 +582,30 @@ void disufq2(double *rsvec, double *isvec,
iv1 = jyi-ixi; iv1 = jyi-ixi;
// iz2 = (n*m-iz1); // iz2 = (n*m-iz1);
// iv2 = (n*m-iv1); // iv2 = (n*m-iv1);
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rsvec[iz1] += tmp1;///rsvec[i][jy+ix] += tmp1; rsvec[iz1] += tmp1;///rsvec[i][jy+ix] += tmp1;
isvec[iz1] += tmp2;///isvec[i][jy+ix] += tmp2; isvec[iz1] += tmp2;///isvec[i][jy+ix] += tmp2;
//rsvec[iz2] += tmp1;///rsvec[i][n*m-(jy+ix)] += tmp1; //rsvec[iz2] += tmp1;///rsvec[i][n*m-(jy+ix)] += tmp1;
//isvec[iz2] -= tmp2;///isvec[i][n*m-(jy-ix)] += tmp2; //isvec[iz2] -= tmp2;///isvec[i][n*m-(jy-ix)] += tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rdvec[iv1] += tmp1; rdvec[iv1] += tmp1;
idvec[iv1] += tmp2; idvec[iv1] += tmp2;
//rdvec[iv2] += tmp1; //rdvec[iv2] += tmp1;
//idvec[iv2] -= tmp2; //idvec[iv2] -= tmp2;
// iv2++; // iv2++;
} }
} }
} }
@ -618,16 +620,16 @@ void disufq2(double *rsvec, double *isvec,
//Visit the MATLAB Central File Exchange for latest version //Visit the MATLAB Central File Exchange for latest version
//http://www.mathworks.com/matlabcentral/fileexchange/3026 //http://www.mathworks.com/matlabcentral/fileexchange/3026
void findrfc3_astm(double *array_ext, double *array_out, int n, int *nout) { void findrfc3_astm(double *array_ext, double *array_out, int n, int *nout) {
double *pr, *po, a[16384], ampl, mean; double *pr, *po, a[16384], ampl, mean;
int tot_num, index, j, cNr1, cNr2; int tot_num, index, j, cNr1, cNr2;
tot_num = n; tot_num = n;
// pointers to the first element of the arrays // pointers to the first element of the arrays
pr = &array_ext[0]; pr = &array_ext[0];
po = &array_out[0]; po = &array_out[0];
// The original rainflow counting by Nieslony, unchanged // The original rainflow counting by Nieslony, unchanged
j = -1; j = -1;
cNr1 = 1; cNr1 = 1;
@ -692,8 +694,8 @@ void
findrfc5_astm(double *array_ext, double *array_t, double *array_out, int n, int *nout) { findrfc5_astm(double *array_ext, double *array_t, double *array_out, int n, int *nout) {
double *pr, *pt, *po, a[16384], t[16384], ampl, mean, period, atime; double *pr, *pt, *po, a[16384], t[16384], ampl, mean, period, atime;
int tot_num, index, j, cNr1, cNr2; int tot_num, index, j, cNr1, cNr2;
// tot_num = mxGetM(array_ext) * mxGetN(array_ext); // tot_num = mxGetM(array_ext) * mxGetN(array_ext);
tot_num = n; tot_num = n;
@ -701,9 +703,9 @@ findrfc5_astm(double *array_ext, double *array_t, double *array_out, int n, int
pr = &array_ext[0]; pr = &array_ext[0];
pt = &array_t[0]; pt = &array_t[0];
po = &array_out[0]; po = &array_out[0];
// array_out = mxCreateDoubleMatrix(5, tot_num-1, mxREAL); // array_out = mxCreateDoubleMatrix(5, tot_num-1, mxREAL);
// The original rainflow counting by Nieslony, unchanged // The original rainflow counting by Nieslony, unchanged
j = -1; j = -1;
cNr1 = 1; cNr1 = 1;

Write Preview
Loading…
Cancel
Save