#ifndef PI
#define PI 3.1415927
#endif

#ifndef G
#define G 9.81
#endif

double jonswap(double omega, double hs, double tp, double gamma)
/*----------------------------------------------------------------------
 *
 * JONSWAP or Pierson-Moskowitz spectrum.
 *
 * For Pierson-Moskowitz, set gamma=1 and the following peak period:
 *    ITTC standard:
 *       The period T1=2*pi*m0/m1 of the spectrum's centroid
 *       is considered input value. The peak period Tp is
 *       Tp = 1.2965 * T1
 *    Bureau Veritas standard:
 *       The mean zero-upcrossing period T0 = 2*pi*sqrt(m0/m2)
 *       is considered input value. The peak period Tp is
 *       Tp = 1.408 * T0
 *
 * Input parameters:
 *    double omega   circular frequency [rad/s]
 *    double hs      significant wave height [m]
 *    double tp      peak period [s]
 *    double gamma   shape factor, must be >= 1; typically 3.3
 *
 * Return value:     energy spectrum [m^2 * s]
 *
 * The Fourier spectrum is:
 *   sqrt(jswap(...)/(2*domega))*cexp(cmplx(0.,2*pi*rnd()))
 *
 * The complex Fourier coefficients are domega * Fourier spectrum, i.e.
 *   sqrt(jswap(...)*domega/2)*cexp(cmplx(0.,2*pi*rnd()))
 *
 * The wave amplitude is twice the amount of the complex Fourier coefficient.
 *
 * Terms of use:
 *    You are not allowed to claim that this code is the work of anybody else
 *    than Henning Weede.
 *    The author is not liable for any damage caused by this code.
 *    The author has no duty to waist his time answering questions related
 *    to this code.
 *    The GNU Public License applies.
 *
 * Author: Henning Weede
 *----------------------------------------------------------------------
 */{

   double a,a0,b,omgp,sigma,tmp,gampol;

   gampol=(gamma-1.)*0.036/2.3;
   if(gampol>0.036) gampol=0.036;

   a0=3.249e-3 * pow(2*PI,4.+2*gampol) *G*G;
   b=1.25*16*PI*PI*PI*PI;
   a=a0*pow(hs/(tp*tp),gampol)*(1.-0.298*log(gamma));
   omgp=2*PI/tp;
   sigma=(omega<omgp)?0.07:0.09;
   tmp=(omega-omgp)/(sigma*omgp);
   tmp=0.5*tmp*tmp;
   if(tmp<30.) a=a*pow(gamma,exp(-tmp));
   tmp=tp*omega;
   tmp=tmp*tmp*tmp*tmp;
   if(tmp>b/30.)
      {
      tmp=b/tmp;
      return(a*hs*hs/(tp*tp*tp*tp*omega*omega*omega*omega*omega)*exp(-tmp));
      }
   else
      return(0.);
/*
   double t1,b,p,to4,s;
   if(omega<0.001) return(0.);
   t1=tp/1.2965;
   b=(omega*t1>=5.32)?0.09:0.07;
   p=omega*t1/5.32-1.; p=p*p/(2*b*b); p=(p<30.)?exp(-p):0.;
   to4=t1*omega; to4=to4*to4*to4*to4;
   s=32.29*PI*PI*PI/to4; s=(s<30.)?exp(-s):0.;
   s=s*pow(gamma,p)*5.32*PI*PI*PI*hs*hs/(to4*omega);
   return(s);
*/
   } /* jonswap() */