Renamed c_codes -> c_library
Per.Andreas.Brodtkorb
10 years ago
parent
2a0efa5ea3
commit
1d006de413
@ -1,98 +0,0 @@
|
|||||||
"""
|
|
||||||
f2py c_library.pyf c_functions.c -c
|
|
||||||
|
|
||||||
See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW
|
|
||||||
"""
|
|
||||||
import os
|
|
||||||
import sys
|
|
||||||
|
|
||||||
def which(program):
|
|
||||||
"""
|
|
||||||
Test if program exists
|
|
||||||
======================
|
|
||||||
|
|
||||||
In order to test if a certain executable exists, it will search for the
|
|
||||||
program name in the environment variables.
|
|
||||||
If program is a full path to an executable, it will check it exists
|
|
||||||
|
|
||||||
Copied from:
|
|
||||||
http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python/
|
|
||||||
It is supposed to mimic the UNIX command "which"
|
|
||||||
"""
|
|
||||||
|
|
||||||
def is_exe(fpath):
|
|
||||||
return os.path.exists(fpath) and os.access(fpath, os.X_OK)
|
|
||||||
|
|
||||||
fpath, unused_fname = os.path.split(program)
|
|
||||||
if fpath:
|
|
||||||
if is_exe(program):
|
|
||||||
return program
|
|
||||||
else:
|
|
||||||
for path in os.environ["PATH"].split(os.pathsep):
|
|
||||||
exe_file = os.path.join(path, program)
|
|
||||||
if is_exe(exe_file):
|
|
||||||
return exe_file
|
|
||||||
|
|
||||||
return None
|
|
||||||
|
|
||||||
def f2py_call_str():
|
|
||||||
# regardless of platform, try to figure out which f2py call is in the path
|
|
||||||
# define possible options
|
|
||||||
|
|
||||||
# on Arch Linux, python and f2py are the calls corresponding to python 3
|
|
||||||
# and python2/f2py2 for python 2
|
|
||||||
# other Linux versions might still use python/f2py for python 2
|
|
||||||
if os.path.basename(sys.executable).endswith('2'):
|
|
||||||
for k in ('f2py2','f2py2.6','f2py2.7',):
|
|
||||||
# if we found the f2py path, no need to look further
|
|
||||||
if which(k):
|
|
||||||
f2py_call = k
|
|
||||||
f2py_path = which(k)
|
|
||||||
break
|
|
||||||
# on Windows and other Linux using python/f2py
|
|
||||||
else:
|
|
||||||
for k in ('f2py','f2py2.6','f2py2.7','f2py.py',):
|
|
||||||
# if we found the f2py path, no need to look further
|
|
||||||
if which(k):
|
|
||||||
f2py_call = k
|
|
||||||
f2py_path = which(k)
|
|
||||||
break
|
|
||||||
|
|
||||||
try:
|
|
||||||
print 'found f2py in:', f2py_path
|
|
||||||
return f2py_call
|
|
||||||
# raise exception if f2py is not found, f2py_path variable will not exist
|
|
||||||
except NameError:
|
|
||||||
raise UserWarning, \
|
|
||||||
'Couldn\'t locate f2py. Should be part of NumPy installation.'
|
|
||||||
# # this might vary among specific cases: f2py, f2py2.7, f2py3.2, ...
|
|
||||||
# # TODO: more robust approach, find out what f2py is in the users path
|
|
||||||
# if os.name == 'posix':
|
|
||||||
# compile_format = 'f2py2.6 %s %s -c'
|
|
||||||
#
|
|
||||||
# # Install microsoft visual c++ .NET 2003 and run the following to build the module:
|
|
||||||
# elif os.name == 'nt':
|
|
||||||
# # compile_format = 'f2py.py %s %s -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71'
|
|
||||||
# compile_format = 'f2py.py %s %s -c'
|
|
||||||
#
|
|
||||||
# # give an Error for other OS-es
|
|
||||||
# else:
|
|
||||||
# raise UserWarning, \
|
|
||||||
# 'Untested platform:', os.name
|
|
||||||
|
|
||||||
def compile_all():
|
|
||||||
f2py_call = f2py_call_str()
|
|
||||||
print '='*75
|
|
||||||
print 'compiling c_codes'
|
|
||||||
print '='*75
|
|
||||||
|
|
||||||
compile_format = f2py_call + ' %s %s -c'
|
|
||||||
|
|
||||||
pyfs = ('c_library.pyf',)
|
|
||||||
files =('c_functions.c',)
|
|
||||||
|
|
||||||
for pyf,file in zip(pyfs,files):
|
|
||||||
os.system(compile_format % (pyf,file))
|
|
||||||
|
|
||||||
if __name__=='__main__':
|
|
||||||
compile_all()
|
|
||||||
@ -1,16 +0,0 @@
|
|||||||
import os
|
|
||||||
|
|
||||||
def compile_all():
|
|
||||||
# Install gfortran and run the following to build the module:
|
|
||||||
#compile_format = 'f2py %s %s -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71'
|
|
||||||
|
|
||||||
# Install microsoft visual c++ .NET 2003 and run the following to build the module:
|
|
||||||
compile_format = 'f2py %s %s -c'
|
|
||||||
pyfs = ('rfc.pyf','diffsumfunq.pyf')
|
|
||||||
files =('findrfc.c','disufq1.c')
|
|
||||||
|
|
||||||
for pyf,file in zip(pyfs,files):
|
|
||||||
os.system(compile_format % (pyf,file))
|
|
||||||
|
|
||||||
if __name__=='__main__':
|
|
||||||
compile_all()
|
|
||||||
Binary file not shown.
@ -1,30 +0,0 @@
|
|||||||
! File diffsumfunq.pyf
|
|
||||||
python module diffsumfunq
|
|
||||||
interface
|
|
||||||
subroutine disufq(rvec, ivec, rA, iA, w, kw, h, g,nmin,nmax, m, n)
|
|
||||||
intent(c) disufq ! disufq is a C function
|
|
||||||
intent(c) ! all disufq arguments are considered as C based
|
|
||||||
!integer intent(hide), depend(rA),check(n*m==len(iA)) :: n=len(rA)/m
|
|
||||||
!integer intent(hide), depend(rA), check(m==shape(iA,1)) :: m=shape(rA,1)
|
|
||||||
double precision dimension(n*m), intent(in) :: rA, iA ! input array
|
|
||||||
double precision dimension(n/2+1), intent(in) :: w, kw ! input array
|
|
||||||
double precision intent(in) :: h, g
|
|
||||||
integer intent(in) :: nmin, nmax
|
|
||||||
double precision dimension(n*m), intent(out) :: rvec, ivec ! output array,
|
|
||||||
end subroutine disufq
|
|
||||||
|
|
||||||
subroutine disufq2(rsvec, isvec,rdvec, idvec, rA, iA, w, kw, h, g,nmin,nmax, m, n)
|
|
||||||
intent(c) disufq2 ! disufq2 is a C function
|
|
||||||
intent(c) ! all disufq2 arguments are considered as C based
|
|
||||||
!integer intent(hide), depend(rA),check(n*m==len(iA)) :: n=len(rA)/m
|
|
||||||
!integer intent(hide), depend(rA), check(m==shape(iA,1)) :: m=shape(rA,1)
|
|
||||||
double precision dimension(n*m), intent(in) :: rA, iA ! input array
|
|
||||||
double precision dimension(n/2+1), intent(in) :: w, kw ! input array
|
|
||||||
double precision intent(in) :: h, g
|
|
||||||
integer intent(in) :: nmin, nmax
|
|
||||||
double precision dimension(n*m), intent(out) :: rsvec, isvec, rdvec, idvec ! output array,
|
|
||||||
|
|
||||||
end subroutine disufq
|
|
||||||
|
|
||||||
end interface
|
|
||||||
end python module diffsumfunq
|
|
||||||
@ -1,30 +0,0 @@
|
|||||||
! File diffsumfunq.pyf
|
|
||||||
python module diffsumfunq
|
|
||||||
interface
|
|
||||||
subroutine disufq(rvec, ivec, rA, iA, w, kw, h, g,nmin,nmax, m, n)
|
|
||||||
intent(c) disufq ! disufq is a C function
|
|
||||||
intent(c) ! all disufq arguments are considered as C based
|
|
||||||
integer intent(hide), depend(rA),check(n==shape(iA,0)) :: n=shape(rA,0)
|
|
||||||
integer intent(hide), depend(rA), check(m==shape(iA,1)) :: m=shape(rA,1)
|
|
||||||
double precision dimension(n,m), intent(in) :: rA, iA ! input array
|
|
||||||
double precision dimension(n), intent(in) :: w, kw ! input array
|
|
||||||
double precision intent(in) :: h, g
|
|
||||||
integer intent(in) :: nmin, nmax
|
|
||||||
double precision dimension(n,m), intent(out) :: rvec, ivec ! output array,
|
|
||||||
end subroutine disufq
|
|
||||||
|
|
||||||
subroutine disufq2(rsvec, isvec,rdvec, idvec, rA, iA, w, kw, h, g,nmin,nmax, m, n)
|
|
||||||
intent(c) disufq2 ! disufq2 is a C function
|
|
||||||
intent(c) ! all disufq2 arguments are considered as C based
|
|
||||||
integer intent(hide), depend(rA),check(n==shape(iA,0)) :: n=shape(rA,0)
|
|
||||||
integer intent(hide), depend(rA), check(m==shape(iA,1)) :: m=shape(rA,1)
|
|
||||||
double precision dimension(n,m), intent(in) :: rA, iA ! input array
|
|
||||||
double precision dimension(n), intent(in) :: w, kw ! input array
|
|
||||||
double precision intent(in) :: h, g
|
|
||||||
integer intent(in) :: nmin, nmax
|
|
||||||
double precision dimension(n,m), intent(out) :: rsvec, isvec, rdvec, idvec ! output array,
|
|
||||||
|
|
||||||
end subroutine disufq
|
|
||||||
|
|
||||||
end interface
|
|
||||||
end python module diffsumfunq
|
|
||||||
@ -1,446 +0,0 @@
|
|||||||
#include "math.h"
|
|
||||||
/*
|
|
||||||
* DISUFQ Is an internal function to spec2nlsdat
|
|
||||||
*
|
|
||||||
* 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).
|
|
||||||
* 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
|
|
||||||
* greater than zero.
|
|
||||||
* 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
|
|
||||||
* 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)).
|
|
||||||
*
|
|
||||||
* Install gfortran and run the following to build the module:
|
|
||||||
* f2py diffsumfunq.pyf disufq1.c -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
|
|
||||||
*
|
|
||||||
* 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 Epij, Edij;
|
|
||||||
double tmp1, tmp2, tmp3, tmp4, kfact;
|
|
||||||
double w1, w2, kw1, kw2, Cg;
|
|
||||||
double rrA, iiA, riA, irA;
|
|
||||||
int i,jy,ix,iz1,iv1,ixi,jyi;
|
|
||||||
//int iz2, iv2;
|
|
||||||
//Initialize rvec and ivec to zero
|
|
||||||
for (ix=0;ix<n*m;ix++) {
|
|
||||||
rvec[ix] = 0.0;
|
|
||||||
ivec[ix] = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
// kfact is set to 2 in order to exploit the symmetry.
|
|
||||||
// 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;
|
|
||||||
if (h>10000){ /* deep water /Inifinite water depth */
|
|
||||||
for (ix = nmin-1;ix<nmax;ix++) {
|
|
||||||
ixi = ix*m;
|
|
||||||
iz1 = 2*ixi;
|
|
||||||
//iz2 = n*m-ixi;
|
|
||||||
kw1 = kw[ix];
|
|
||||||
Epij = kw1;
|
|
||||||
for (i=0;i<m;i++,ixi++,iz1++) {
|
|
||||||
rrA = rA[ixi]*rA[ixi]; ///
|
|
||||||
iiA = iA[ixi]*iA[ixi]; ///
|
|
||||||
riA = rA[ixi]*iA[ixi]; ///
|
|
||||||
|
|
||||||
/// Sum frequency effects along the diagonal
|
|
||||||
tmp1 = kfact*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*riA*Epij;
|
|
||||||
rvec[iz1] += tmp1;
|
|
||||||
ivec[iz1] += tmp2;
|
|
||||||
|
|
||||||
//rvec[iz2] += tmp1;
|
|
||||||
//ivec[iz2] -= tmp2;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/// Difference frequency effects are zero along the diagonal
|
|
||||||
/// and are thus not contributing to the mean.
|
|
||||||
}
|
|
||||||
for (jy = ix+1;jy<nmax;jy++){
|
|
||||||
kw2 = kw[jy];
|
|
||||||
Epij = 0.5*(kw2 + kw1);
|
|
||||||
Edij = -0.5*(kw2 - kw1);
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
jyi = jy*m;
|
|
||||||
iz1 = ixi+jyi;
|
|
||||||
iv1 = jyi-ixi;
|
|
||||||
//iz2 = (n*m-iz1);
|
|
||||||
//iv2 = (n*m-iv1);
|
|
||||||
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
|
|
||||||
|
|
||||||
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[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];
|
|
||||||
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
|
|
||||||
|
|
||||||
/* Sum frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*(riA+irA)*Epij;
|
|
||||||
rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
|
|
||||||
ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
|
|
||||||
//rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
|
|
||||||
//ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2;
|
|
||||||
// iz2++;
|
|
||||||
|
|
||||||
/* Difference frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA+iiA)*Edij;
|
|
||||||
tmp2 = kfact*2.0*(riA-irA)*Edij;
|
|
||||||
|
|
||||||
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
|
|
||||||
ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
|
|
||||||
|
|
||||||
//rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
|
|
||||||
//ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
|
|
||||||
//iv2++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else{ /* Finite water depth */
|
|
||||||
for (ix = nmin-1;ix<nmax;ix++) {
|
|
||||||
kw1 = kw[ix];
|
|
||||||
w1 = w[ix];
|
|
||||||
tmp1 = tanh(kw1*h);
|
|
||||||
/// Cg, wave group velocity
|
|
||||||
Cg = 0.5*g*(tmp1 + kw1*h*(1.0- tmp1*tmp1))/w1; /// OK
|
|
||||||
tmp1 = 0.5*g*(kw1/w1)*(kw1/w1);
|
|
||||||
tmp2 = 0.5*w1*w1/g;
|
|
||||||
tmp3 = g*kw1/(w1*Cg);
|
|
||||||
|
|
||||||
if (kw1*h<300.0){
|
|
||||||
tmp4 = kw1/sinh(2.0*kw1*h);
|
|
||||||
}
|
|
||||||
else{ // To ensure sinh does not overflow.
|
|
||||||
tmp4 = 0.0;
|
|
||||||
}
|
|
||||||
// Difference frequency effects finite water depth
|
|
||||||
Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK
|
|
||||||
|
|
||||||
// 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
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
iz1 = 2*ixi;
|
|
||||||
//iz2 = n*m-ixi;
|
|
||||||
for (i=0;i<m;i++,ixi++,iz1++) {
|
|
||||||
|
|
||||||
rrA = rA[ixi]*rA[ixi]; ///
|
|
||||||
iiA = iA[ixi]*iA[ixi]; ///
|
|
||||||
riA = rA[ixi]*iA[ixi]; ///
|
|
||||||
|
|
||||||
|
|
||||||
/// Sum frequency effects along the diagonal
|
|
||||||
rvec[iz1] += kfact*(rrA-iiA)*Epij;
|
|
||||||
ivec[iz1] += kfact*2.0*riA*Epij;
|
|
||||||
//rvec[iz2] += kfact*(rrA-iiA)*Epij;
|
|
||||||
//ivec[iz2] -= kfact*2.0*riA*Epij;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/// Difference frequency effects along the diagonal
|
|
||||||
/// are only contributing to the mean
|
|
||||||
rvec[i] += 2.0*(rrA+iiA)*Edij;
|
|
||||||
}
|
|
||||||
for (jy = ix+1;jy<nmax;jy++) {
|
|
||||||
// w1 = w[ix];
|
|
||||||
// kw1 = kw[ix];
|
|
||||||
w2 = w[jy];
|
|
||||||
kw2 = kw[jy];
|
|
||||||
tmp1 = g*(kw1/w1)*(kw2/w2);
|
|
||||||
tmp2 = 0.5/g*(w1*w1+w2*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));
|
|
||||||
Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* 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));
|
|
||||||
tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h));
|
|
||||||
Edij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
jyi = jy*m;
|
|
||||||
iz1 = ixi+jyi;
|
|
||||||
iv1 = jyi-ixi;
|
|
||||||
// iz2 = (n*m-iz1);
|
|
||||||
// iv2 = n*m-iv1;
|
|
||||||
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
|
|
||||||
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[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];
|
|
||||||
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
|
|
||||||
|
|
||||||
/* Sum frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*(riA+irA)*Epij;
|
|
||||||
rvec[iz1] += tmp1;///rvec[i][jy+ix] += tmp1;
|
|
||||||
ivec[iz1] += tmp2;///ivec[i][jy+ix] += tmp2;
|
|
||||||
//rvec[iz2] += tmp1;///rvec[i][n*m-(jy+ix)] += tmp1;
|
|
||||||
//ivec[iz2] -= tmp2;///ivec[i][n*m-(jy+ix)] -= tmp2;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/* Difference frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA+iiA)*Edij;
|
|
||||||
tmp2 = kfact*2.0*(riA-irA)*Edij;
|
|
||||||
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
|
|
||||||
ivec[iv1] += tmp2;///ivec[i][jy-ix] -= tmp2;
|
|
||||||
|
|
||||||
//rvec[iv2] += tmp1;
|
|
||||||
//ivec[iv2] -= tmp2;
|
|
||||||
//iv2++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
//return i;
|
|
||||||
}
|
|
||||||
/*
|
|
||||||
* DISUFQ2 Is an internal function to spec2nlsdat
|
|
||||||
*
|
|
||||||
* CALL: disufq2(rsvec,isvec,rdvec,idvec,rA,iA, w,kw,h,g,nmin,nmax,m,n)
|
|
||||||
*
|
|
||||||
* rsvec, isvec = real and imaginary parts of the sum frequency
|
|
||||||
* effects (size m X n).
|
|
||||||
* rdvec, idvec = real and imaginary parts of the difference frequency
|
|
||||||
* effects (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
|
|
||||||
* greater than zero.
|
|
||||||
* 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
|
|
||||||
* 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
|
|
||||||
* a simple 1D Fourier transform, real(FFT(svec+dvec)).
|
|
||||||
*
|
|
||||||
*
|
|
||||||
* This is a MEX-file for MATLAB.
|
|
||||||
* by Per Andreas Brodtkorb 15.08.2001
|
|
||||||
* 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 Epij, Edij;
|
|
||||||
double tmp1, tmp2, tmp3, tmp4, kfact;
|
|
||||||
double w1, w2, kw1, kw2, Cg;
|
|
||||||
double rrA, iiA, riA, irA;
|
|
||||||
int i,jy,ix,iz1,iv1,ixi,jyi;
|
|
||||||
//int iz2,iv2
|
|
||||||
|
|
||||||
//Initialize rvec and ivec to zero
|
|
||||||
for (ix=0;ix<n*m;ix++) {
|
|
||||||
rsvec[ix] = 0.0;
|
|
||||||
isvec[ix] = 0.0;
|
|
||||||
rdvec[ix] = 0.0;
|
|
||||||
idvec[ix] = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
// kfact is set to 2 in order to exploit the symmetry.
|
|
||||||
// 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;
|
|
||||||
if (h>10000){ /* deep water /Inifinite water depth */
|
|
||||||
for (ix = nmin-1;ix<nmax;ix++) {
|
|
||||||
ixi = ix*m;
|
|
||||||
iz1 = 2*ixi;
|
|
||||||
//iz2 = n*m-ixi;
|
|
||||||
kw1 = kw[ix];
|
|
||||||
Epij = kw1;
|
|
||||||
for (i=0;i<m;i++,ixi++,iz1++) {
|
|
||||||
rrA = rA[ixi]*rA[ixi]; ///
|
|
||||||
iiA = iA[ixi]*iA[ixi]; ///
|
|
||||||
riA = rA[ixi]*iA[ixi]; ///
|
|
||||||
|
|
||||||
/// Sum frequency effects along the diagonal
|
|
||||||
tmp1 = kfact*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*riA*Epij;
|
|
||||||
rsvec[iz1] += tmp1;
|
|
||||||
isvec[iz1] += tmp2;
|
|
||||||
|
|
||||||
//rsvec[iz2] += tmp1;
|
|
||||||
//isvec[iz2] -= tmp2;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/// Difference frequency effects are zero along the diagonal
|
|
||||||
/// and are thus not contributing to the mean.
|
|
||||||
}
|
|
||||||
for (jy = ix+1;jy<nmax;jy++){
|
|
||||||
kw2 = kw[jy];
|
|
||||||
Epij = 0.5*(kw2 + kw1);
|
|
||||||
Edij = -0.5*(kw2 - kw1);
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
jyi = jy*m;
|
|
||||||
iz1 = ixi+jyi;
|
|
||||||
iv1 = jyi-ixi;
|
|
||||||
//iz2 = (n*m-iz1);
|
|
||||||
//iv2 = (n*m-iv1);
|
|
||||||
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
|
|
||||||
|
|
||||||
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[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];
|
|
||||||
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
|
|
||||||
|
|
||||||
/* Sum frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*(riA+irA)*Epij;
|
|
||||||
rsvec[iz1] += tmp1; ///rvec[i][ix+jy] += tmp1;
|
|
||||||
isvec[iz1] += tmp2; ///ivec[i][ix+jy] += tmp2;
|
|
||||||
//rsvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
|
|
||||||
//isvec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] += tmp2;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/* Difference frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA+iiA)*Edij;
|
|
||||||
tmp2 = kfact*2.0*(riA-irA)*Edij;
|
|
||||||
|
|
||||||
rdvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
|
|
||||||
idvec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
|
|
||||||
|
|
||||||
//rdvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
|
|
||||||
//idvec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
|
|
||||||
// iv2++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else{ /* Finite water depth */
|
|
||||||
for (ix = nmin-1;ix<nmax;ix++) {
|
|
||||||
kw1 = kw[ix];
|
|
||||||
w1 = w[ix];
|
|
||||||
tmp1 = tanh(kw1*h);
|
|
||||||
/// Cg, wave group velocity
|
|
||||||
Cg = 0.5*g*(tmp1 + kw1*h*(1.0- tmp1*tmp1))/w1; /// OK
|
|
||||||
tmp1 = 0.5*g*(kw1/w1)*(kw1/w1);
|
|
||||||
tmp2 = 0.5*w1*w1/g;
|
|
||||||
tmp3 = g*kw1/(w1*Cg);
|
|
||||||
|
|
||||||
if (kw1*h<300.0){
|
|
||||||
tmp4 = kw1/sinh(2.0*kw1*h);
|
|
||||||
}
|
|
||||||
else{ // To ensure sinh does not overflow.
|
|
||||||
tmp4 = 0.0;
|
|
||||||
}
|
|
||||||
// Difference frequency effects finite water depth
|
|
||||||
Edij = (tmp1-tmp2+tmp3)/(1.0-g*h/(Cg*Cg))-tmp4; /// OK
|
|
||||||
|
|
||||||
// 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
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
iz1 = 2*ixi;
|
|
||||||
//iz2 = n*m-ixi;
|
|
||||||
for (i=0;i<m;i++,ixi++,iz1++) {
|
|
||||||
|
|
||||||
rrA = rA[ixi]*rA[ixi]; ///
|
|
||||||
iiA = iA[ixi]*iA[ixi]; ///
|
|
||||||
riA = rA[ixi]*iA[ixi]; ///
|
|
||||||
|
|
||||||
|
|
||||||
/// Sum frequency effects along the diagonal
|
|
||||||
rsvec[iz1] += kfact*(rrA-iiA)*Epij;
|
|
||||||
isvec[iz1] += kfact*2.0*riA*Epij;
|
|
||||||
//rsvec[iz2] += kfact*(rrA-iiA)*Epij;
|
|
||||||
//isvec[iz2] -= kfact*2.0*riA*Epij;
|
|
||||||
|
|
||||||
/// Difference frequency effects along the diagonal
|
|
||||||
/// are only contributing to the mean
|
|
||||||
//printf(" %f \n",2.0*(rrA+iiA)*Edij);
|
|
||||||
rdvec[i] += 2.0*(rrA+iiA)*Edij;
|
|
||||||
}
|
|
||||||
for (jy = ix+1;jy<nmax;jy++) {
|
|
||||||
// w1 = w[ix];
|
|
||||||
// kw1 = kw[ix];
|
|
||||||
w2 = w[jy];
|
|
||||||
kw2 = kw[jy];
|
|
||||||
tmp1 = g*(kw1/w1)*(kw2/w2);
|
|
||||||
tmp2 = 0.5/g*(w1*w1+w2*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));
|
|
||||||
Epij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* 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));
|
|
||||||
tmp4 = (1.0-g*(kw1-kw2)/(w1-w2)/(w1-w2)*tanh((kw1-kw2)*h));
|
|
||||||
Edij = (tmp1-tmp2+tmp3)/tmp4+tmp2-0.5*tmp1; /* OK */
|
|
||||||
//printf("Edij = %f Epij = %f \n", Edij,Epij);
|
|
||||||
|
|
||||||
ixi = ix*m;
|
|
||||||
jyi = jy*m;
|
|
||||||
iz1 = ixi+jyi;
|
|
||||||
iv1 = jyi-ixi;
|
|
||||||
// iz2 = (n*m-iz1);
|
|
||||||
// iv2 = (n*m-iv1);
|
|
||||||
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
|
|
||||||
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[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];
|
|
||||||
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
|
|
||||||
|
|
||||||
/* Sum frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA-iiA)*Epij;
|
|
||||||
tmp2 = kfact*2.0*(riA+irA)*Epij;
|
|
||||||
rsvec[iz1] += tmp1;///rsvec[i][jy+ix] += tmp1;
|
|
||||||
isvec[iz1] += tmp2;///isvec[i][jy+ix] += tmp2;
|
|
||||||
//rsvec[iz2] += tmp1;///rsvec[i][n*m-(jy+ix)] += tmp1;
|
|
||||||
//isvec[iz2] -= tmp2;///isvec[i][n*m-(jy-ix)] += tmp2;
|
|
||||||
//iz2++;
|
|
||||||
|
|
||||||
/* Difference frequency effects */
|
|
||||||
tmp1 = kfact*2.0*(rrA+iiA)*Edij;
|
|
||||||
tmp2 = kfact*2.0*(riA-irA)*Edij;
|
|
||||||
rdvec[iv1] += tmp1;
|
|
||||||
idvec[iv1] += tmp2;
|
|
||||||
|
|
||||||
//rdvec[iv2] += tmp1;
|
|
||||||
//idvec[iv2] -= tmp2;
|
|
||||||
// iv2++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// return i;
|
|
||||||
}
|
|
||||||
@ -1,53 +0,0 @@
|
|||||||
|
|
||||||
|
|
||||||
/*
|
|
||||||
* 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, double *ind, int n, int info)
|
|
||||||
{ int i,start, ix=0,dcross=0;
|
|
||||||
|
|
||||||
if ( *(y +0)< v){
|
|
||||||
dcross=-1; /* first is a up-crossing*/
|
|
||||||
}
|
|
||||||
if ( *(y +0)> v){
|
|
||||||
dcross=1; /* first is a down-crossing*/
|
|
||||||
}
|
|
||||||
start=0;
|
|
||||||
if ( *(y +0)== v){
|
|
||||||
/* Find out what type of crossing we have next time.. */
|
|
||||||
for (i=1; i<n; i++) {
|
|
||||||
start=i;
|
|
||||||
if ( *(y +i)< v){
|
|
||||||
*(ind + ix) = i; /* first crossing is a down crossing*/
|
|
||||||
ix++;
|
|
||||||
dcross=-1; /* The next crossing is a up-crossing*/
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
if ( *(y +i)> v){
|
|
||||||
*(ind + ix) = i; /* first crossing is a up-crossing*/
|
|
||||||
ix++;
|
|
||||||
dcross=1; /*The next crossing is a down-crossing*/
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (i=start; i<n-1; i++) {
|
|
||||||
if (( (dcross==-1) && (*(y +i)<=h) && (*(y+i+1) > h) ) || ((dcross==1 ) && (*(y +i)>=h) && (*(y+i+1) < h) ) ) {
|
|
||||||
|
|
||||||
*(ind + ix) = i+1 ;
|
|
||||||
ix++;
|
|
||||||
dcross=-dcross;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
info = ix
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,118 +0,0 @@
|
|||||||
#include "math.h"
|
|
||||||
/*
|
|
||||||
* findrfc.c -
|
|
||||||
*
|
|
||||||
* Returns indices to RFC turningpoints of a vector
|
|
||||||
* of turningpoints
|
|
||||||
*
|
|
||||||
* Install gfortran and run the following to build the module:
|
|
||||||
* f2py rfc.pyf findrfc.c -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
|
|
||||||
*
|
|
||||||
* 1998 by Per Andreas Brodtkorb.
|
|
||||||
*/
|
|
||||||
|
|
||||||
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,iy;
|
|
||||||
|
|
||||||
if (*(y1+0)> *(y1+1)){
|
|
||||||
/* if first is a max , ignore the first max*/
|
|
||||||
y=&(*(y1+1));
|
|
||||||
NC=floor((n-1)/2);
|
|
||||||
Tstart=2;
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
y=y1;
|
|
||||||
NC=floor(n/2);
|
|
||||||
Tstart=1;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (NC<1){
|
|
||||||
info = 0;
|
|
||||||
return; /* No RFC cycles*/
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){
|
|
||||||
info = -1;
|
|
||||||
return; /*This is not a sequence of turningpoints, exit */
|
|
||||||
}
|
|
||||||
if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){
|
|
||||||
info=-1;
|
|
||||||
return; /*This is not a sequence of turningpoints, exit */
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
for (i=0; i<NC; i++) {
|
|
||||||
|
|
||||||
Tmi=Tstart+2*i;
|
|
||||||
Tpl=Tstart+2*i+2;
|
|
||||||
xminus=*(y+2*i);
|
|
||||||
xplus=*(y+2*i+2);
|
|
||||||
|
|
||||||
if(i!=0){
|
|
||||||
j=i-1;
|
|
||||||
while((j>=0) && (*(y+2*j+1)<=*(y+2*i+1))){
|
|
||||||
if( (*(y+2*j)<xminus) ){
|
|
||||||
xminus=*(y+2*j);
|
|
||||||
Tmi=Tstart+2*j;
|
|
||||||
} /*if */
|
|
||||||
j--;
|
|
||||||
} /*while j*/
|
|
||||||
} /*if i */
|
|
||||||
if ( xminus >= xplus){
|
|
||||||
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<NC) ) {
|
|
||||||
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++;
|
|
||||||
|
|
||||||
} /*if*/
|
|
||||||
goto L180;
|
|
||||||
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;*/
|
|
||||||
}
|
|
||||||
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 = ix;
|
|
||||||
return ;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
Binary file not shown.
@ -1,14 +0,0 @@
|
|||||||
! File rfc.pyf
|
|
||||||
python module rfc
|
|
||||||
interface
|
|
||||||
subroutine findrfc(y1,hmin, ind, n,info)
|
|
||||||
intent(c) findrfc ! findrfc is a C function
|
|
||||||
intent(c) ! all findrfc arguments are considered as C based
|
|
||||||
integer intent(hide), depend(y1) :: n=len(y1)
|
|
||||||
double precision dimension(n), intent(in) :: y1 ! input array
|
|
||||||
double precision intent(in) :: hmin
|
|
||||||
integer dimension(n), intent(out) :: ind ! output array,
|
|
||||||
integer intent(out) :: info
|
|
||||||
end subroutine findrfc
|
|
||||||
end interface
|
|
||||||
end python module rfc
|
|
||||||
@ -0,0 +1,25 @@
|
|||||||
|
"""
|
||||||
|
f2py c_library.pyf c_functions.c -c
|
||||||
|
|
||||||
|
See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW
|
||||||
|
"""
|
||||||
|
import os
|
||||||
|
import sys
|
||||||
|
from wafo.f2py_tools import f2py_call_str
|
||||||
|
|
||||||
|
def compile_all():
|
||||||
|
f2py_call = f2py_call_str()
|
||||||
|
print '=' * 75
|
||||||
|
print 'compiling c_codes'
|
||||||
|
print '=' * 75
|
||||||
|
|
||||||
|
compile_format = f2py_call + ' %s %s -c'
|
||||||
|
|
||||||
|
pyfs = ('c_library.pyf',)
|
||||||
|
files = ('c_functions.c',)
|
||||||
|
|
||||||
|
for pyf, file_ in zip(pyfs, files):
|
||||||
|
os.system(compile_format % (pyf, file_))
|
||||||
|
|
||||||
|
if __name__ == '__main__':
|
||||||
|
compile_all()
|
||||||
Binary file not shown.
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
Binary file not shown.
Loading…
Reference in New Issue