From 1863b2233fcb540bf07e988a82e24bd6bcade118 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Sat, 13 Feb 2016 22:01:36 +0100 Subject: [PATCH] Restored c_functions.c --- wafo/source/c_library/c_functions.c | 624 ++++++++++++++-------------- 1 file changed, 313 insertions(+), 311 deletions(-) diff --git a/wafo/source/c_library/c_functions.c b/wafo/source/c_library/c_functions.c index d887b20..a9bbdd5 100644 --- a/wafo/source/c_library/c_functions.c +++ b/wafo/source/c_library/c_functions.c @@ -1,11 +1,11 @@ -#include "math.h" +#include "math.h" /* * 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 */ - + /* - * findrfc.c - + * findrfc.c - * * Returns indices to RFC turningpoints of a vector * of turningpoints @@ -15,11 +15,11 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) { double xminus,xplus,Tpl,Tmi,*y,Tstart; - int i,j,ix=0,NC; + int i,j,ix=0,NC,iy; info[0] = 0; - if (*(y1+0)> *(y1+1)){ + if (*(y1+0)> *(y1+1)){ /* if first is a max , ignore the first max*/ - y=&(*(y1+1)); + y=&(*(y1+1)); NC=floor((n-1)/2); Tstart=1; } @@ -28,140 +28,142 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) { NC=floor(n/2); Tstart=0; } - + if (NC<1){ return; /* No RFC cycles*/ } - + if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){ - info[0] = -1; + info[0] = -1; return; /*This is not a sequence of turningpoints, exit */ } if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){ info[0]=-1; return; /*This is not a sequence of turningpoints, exit */ } - - + + for (i=0; i=0) && (*(y+2*j+1)<=*(y+2*i+1))){ - if( (*(y+2*j)=0) && (*(y+2*j+1)<=*(y+2*i+1))){ + if( (*(y+2*j)= xplus){ - if ( (*(y+2*i+1)-xminus) >= hmin){ - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; - } /*if*/ - goto L180; - } - + if ( (*(y+2*i+1)-xminus) >= hmin){ + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ + goto L180; + } + j=i+1; while((j= *(y+2*i+1)) goto L170; - if( (*(y+2*j+2) <= xplus) ){ - xplus=*(y+2*j+2); - Tpl=(Tstart+2*j+2); - }/*if*/ - j++; + if (*(y+2*j+1) >= *(y+2*i+1)) goto L170; + if( (*(y+2*j+2) <= xplus) ){ + xplus=*(y+2*j+2); + Tpl=(Tstart+2*j+2); + }/*if*/ + j++; } /*while*/ - + if ( (*(y+2*i+1)-xminus) >= hmin) { - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ goto L180; - L170: + L170: if (xplus <= xminus ) { - if ( (*(y+2*i+1)-xminus) >= hmin){ - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; - } /*if*/ - /*goto L180;*/ + if ( (*(y+2*i+1)-xminus) >= hmin){ + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ + /*goto L180;*/ } - else{ - if ( (*(y+2*i+1)-xplus) >= hmin) { - *(ind+ix)=(Tstart+2*i+1); - ix++; - *(ind+ix)=Tpl; - ix++; - } /*if*/ + else{ + if ( (*(y+2*i+1)-xplus) >= hmin) { + *(ind+ix)=(Tstart+2*i+1); + ix++; + *(ind+ix)=Tpl; + ix++; + } /*if*/ } /*elseif*/ L180: + iy=i; } /* for i */ info[0] = ix; - return ; + return ; } /* - * findcross.c - + * findcross.c - * * Returns indices to level v crossings of argument vector * * 1998 by Per Andreas Brodtkorb. last modified 23.06-98 - */ + */ void findcross(double *y, double v, int *ind, int n, int *info) { int i,start, ix=0,dcross=0; - start=0; - if ( y[0]< v){ - dcross=-1; /* first is a up-crossing*/ - } - else if ( y[0]> v){ - dcross=1; /* first is a down-crossing*/ - } - else if ( y[0]== v){ - /* Find out what type of crossing we have next time.. */ - for (i=1; i v){ - ind[ix] = i-1; /* first crossing is a up-crossing*/ - ix++; - dcross=1; /*The next crossing is a down-crossing*/ - goto L120; - } - } - } - L120: - for (i=start; i v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { - - ind[ix] = i; - ix++; - dcross=-dcross; - } - } - info[0] = ix; - return; + start=0; + if ( y[0]< v){ + dcross=-1; /* first is a up-crossing*/ + } + else if ( y[0]> v){ + dcross=1; /* first is a down-crossing*/ + } + else if ( y[0]== v){ + /* Find out what type of crossing we have next time.. */ + for (i=1; i v){ + ind[ix] = i-1; /* first crossing is a up-crossing*/ + ix++; + dcross=1; /*The next crossing is a down-crossing*/ + goto L120; + } + } + } + L120: + for (i=start; i v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { + + ind[ix] = i; + ix++; + dcross=-dcross; + } + } + info[0] = ix; + 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) * * 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) * kw = vector with wavenumbers (kw>=0) * h = water depth (h >=0) * 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. - * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is + * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * greater than zero. * m = size(rA,1),size(iA,1) * 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. * The 2'nd order contribution to the Stokes wave is then calculated by * 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 * revised pab 14.03.2002, 01.05.2002 22.07.2002, oct 2008 - */ + */ void disufq(double *rvec, double *ivec, - double *rA, double *iA, - double *w, double *kw, - double h, double g, - int nmin, int nmax, - int m, int n) -{ + double *rA, double *iA, + double *w, double *kw, + double h, double g, + int nmin, int nmax, + int m, int n) +{ double Epij, Edij; double tmp1, tmp2, tmp3, tmp4, kfact; 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. - // 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]. kfact = 2.0; @@ -225,70 +227,70 @@ void disufq(double *rvec, double *ivec, iz1 = 2*ixi; //iz2 = n*m-ixi; kw1 = kw[ix]; - Epij = kw1; - for (i=0;i=0) * kw = vector with wavenumbers (kw>=0) * h = water depth (h >=0) * 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. - * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is + * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * greater than zero. * m = size(rA,1),size(iA,1) * 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 * dvec = rdvec +sqrt(-1)*idvec. * 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. * 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, - double *rdvec, double *idvec, - double *rA, double *iA, - double *w, double *kw, - double h, double g, - int nmin, int nmax, - int m, int n) -{ + double *rdvec, double *idvec, + double *rA, double *iA, + double *w, double *kw, + double h, double g, + int nmin, int nmax, + int m, int n) +{ double Epij, Edij; double tmp1, tmp2, tmp3, tmp4, kfact; 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. - // 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]. kfact = 2.0; @@ -451,70 +453,70 @@ void disufq2(double *rsvec, double *isvec, iz1 = 2*ixi; //iz2 = n*m-ixi; kw1 = kw[ix]; - Epij = kw1; - for (i=0;i