// Persistence of Vision Ray Tracer Include File
// File: waves.inc
// Desc: water wave simulation      
//       based on technique described in
//       http://www.first.gmd.de/persons/bwalter/html/report/report.html
//       and implementation from issue 12/2001 of 'PC-Magazin'
// Date: November 2001
// Auth: Christoph Hormann <chris_hormann@gmx.de>
// 
// Version 0.1 (experimantal)

#include "gamma.inc"  

#declare gravity = 9.81;
#declare MAX_WAVES = 40;
#declare xalpha = 0.008;  

#declare Seed=seed(9);

#macro directionEnergy(F, Theta)
  #local Temp = cos(Theta*0.5);
  #local Temp = Temp^2;
  
  #if (F >= freqPeak)
    #local s_sm = (F/freqPeak)^-2.5;
    #local sm = 9.77;
  #else
    #local s_sm = (F/freqPeak)^5.0;
    #local sm = 6.97;
  #end 

  #local S = s_sm * sm;   
  (gamma( S+1 ) * Temp^S / (2.0*sqrt(pi)*gamma(S+0.5)) )
#end 

#macro energylDFreq(F)
  ( xalpha * gravity^2 / ( (2.0*pi)^4 * F^5 ) * exp((-5.0/4.0)* (F/freqPeak)^-4.0 ) )
#end 

#macro energy2DFreqAngle(F, Theta)   
  ( energylDFreq( F ) * directionEnergy( F, Theta ) )
#end

#macro waveAmplitude( F, Theta, K )  
 #local Omega = 10.0;
 ( sqrt( energy2DFreqAngle( F, Theta ) * gravity * pi^2 / ( K * F * Omega ) ) )
#end



#macro Waves(windDirection, windSpeed, time, Scale)

  #local freqPeak = 0.13 * gravity / windSpeed;

  ( 
    
    #debug "generating function...\n"

    #local nWaves = 0;
          
    #while ( nWaves < MAX_WAVES )

      #local Frequency = freqPeak + rand(Seed) - 0.5;

      #if ( Frequency > 0.0 )
        
        #local Wave_lambda = 2*pi / ( (2.0 * pi * Frequency)^2 / gravity );
        #local Wave_k = 2*pi / Wave_lambda;
        #local Wave_omega = sqrt(gravity * Wave_k );
        #local Wave_freq = Wave_omega / ( 2*pi );
        #local Wave_periode = 1.0 / Wave_freq;
                                                                              
        #local Wave_direction = ( rand(Seed)*2 ) * pi;
        #local Wave_dirX = cos( Wave_direction ) * 0.5;
        #local Wave_dirY = sin( Wave_direction ) * 0.5;

        #local phi0 = rand(Seed)*2*pi;
        #local a0 = waveAmplitude( Wave_freq, Wave_direction - windDirection, Wave_k );
        #local Wave_amplitude = a0 * cos( phi0 );
        #local Wave_phase = a0 * sin( phi0 );

        #if ( abs(a0) >= 0.0001 )
          #local nWaves = nWaves + 1;
        #end  
        
        // --- Function code ---
        
        (
        Wave_amplitude *            
        sin( Wave_phase + Wave_omega*Wave_freq*( time + Wave_dirX*x*Scale + Wave_dirY*y*Scale ) )
        )
        
        #if ( nWaves < MAX_WAVES )         
          +        
        #end
            
      #end
    #end 
  
  )
  
#end