#version 3.6;

#include "colors.inc"
#include "functions.inc"

global_settings {
  // Not using "assumed_gamma". This way, output of wake_waves function is
  // lineairly proportional to the resulting grey-scaled colors in the output
  // image.
  // assumed_gamma 1.0
}

// ----------------------------------------

#local Zoom = 10;
camera {
	orthographic 
	location < 0, 0, -5>
	right x * Zoom
	up y * Zoom
	look_at < 0, 0.0, 0.0>
}


#macro PlotFunction(Fun, MinX, MaxX, Inc, Size, Color)
  #local X = MinX;
  #while (X < MaxX)
    sphere {
      <X, Fun(X), 0>, Size
      texture {
        pigment { color Color }
        finish { ambient 1 }
      }
    }
    #local X = X + Inc;
  #end
#end


#local Chebyshev_dist = function(x, y, x2, y2) {
	   max(abs(x-x2), abs(y-y2))
  }

#local pos_cos = function(x) {
    0.5 + 0.5 * cos(x)
  }

// Function "x modulus a" with modulus defined such that it is always
// positive.
#local abs_mod = function(x, a) {
    select(x, a-mod(-x, a), mod(x, a))
  }

#local delta = function(x) {
    select(abs(x)-1, 1-abs(x), 0)
  }

#local even = function(x) {
    select( abs_mod(x, 2) - 1, 0, 1 )
  }
#local pos_even = function(x) {
    select( x, 0, select( mod(x, 2) - 1, 0, 1 ) )
  }

#local checker_ = function(x, y, a) {
    even( even(x*a) +  even(y * a))
  }

// Symmetrical dampening function. Unit scaling at x=0, approaching zero for
// large values of x. Fall-off increases for larger values of a. 
#local damp_sym = function(x, a) {
    pos_cos(2 * atan2(a * x, 1))
  }

// A-symmetrical dampening function. Fall-off for negative x controlled
// by a, for positive x controlled by b.
#local damp_asym = function(x, a, b) {
    select(x, damp_sym(x, a), damp_sym(x, b))
  }

// Dampening of possitive x, fall-off controlled by b. For negative x, 
// zero, except for -a < x < 0, where it gradually goes from zero to
// one.
#local damp_half = function(x, a, b) {
    select(x, 
      select(x+a, 0, pos_cos(pi * x/a)),
      damp_sym(x, b))
  }

// Wave
#local wave = function(x) {
    pos_cos(x) * atan2(x, 1)
  }

// Pre: x in [-b, b] range.
// Output increasing from -b to b, but not lineair, but polynomial raised
// to the power a. 
#local scale_ = function(x, a, b) {
    select(x,
      b * (pow((x+b)/b, a) - 1), 
      b * (1 - pow((b-x)/b, a)) )
  }

// Repeating waves (period 2*pi) with a nearly flat bit inbetween, the
// size of which is controlled by a.
#local waves_ = function(x, a) {
    -wave( scale_( abs_mod(x + pi, pi * 2) - pi, a, pi ) )
  }

// Repeating waves (period 2*pi), simply sinus-shaped
#local waves2_ = function(x) {
    sin(x)
  }

// The distance of the wake's peak (if the wake front is extended in
// front of the ship) from the original of the waves' coordinate system.
//
// w: width of (this half of) the wake at (0, 0) in the wake's coordinate
// system.
// wake_angle: angle of the wake (in radians) 
#local peak_dist = function(w, wake_angle) {
    w / sin( wake_angle )
  }

#local phase_shift = 0;

#local phase_outside = function(x, y, w, wake_angle) {
    x + 1.5 * atan2(y, 1) + 1.5 * y + peak_dist(w, wake_angle) 
      + phase_shift - w / tan(wake_angle)
  }

#local phase_outside2 = function(x, y, w, wake_angle) {
    phase_outside(x, y / (1 + 0 * 0.1 * x), w, wake_angle )
  } 

// The alpha coordinate (in radians) of the inner-wave angular coordinate 
// system given the coordinate in the wake's cartesian coordinate system.
#local alpha_ = function(x, y, w, wake_angle) {
    atan2(y, x + peak_dist(w, wake_angle))
  }

// The radius coordinate of the inner-wave angular coordinate system
// given the coordinate in the wake's cartesian coordinate system.
#local radius_ = function(x, y, w, wake_angle) {
    sqrt( pow(x + peak_dist(w, wake_angle), 2) 
          + pow(y, 2) )
  }

#local phase_inside = function(x, y, w, wake_angle, a, b) {
    (phase_shift + radius_(x, y, w, wake_angle) - w / tan(wake_angle) )
    * (1 + pow( a, -b)
         - pow( (wake_angle - alpha_(x, y, w, wake_angle)) / (a * wake_angle)
                , b)) 
  }

#local wave_freq = 3;
#local phase_ = function(x, y, w, wake_angle) {
    -pi * 0 + wave_freq * pi * select(y, 
                 phase_outside2(x, y, w, wake_angle), 
                 phase_inside(x, y, w, wake_angle, 1.5, 2.9 ))
  }

// Averages the phase_ function around y=0 to get a more gradual transition. 
// Parameter a controls the distance of the two points that are averaged,
// and b controls the range of y-values within which averaging takes place.
#local phase2_ = function(x, y, w, wake_angle, a, b) {
    ( phase_(x, y - a * delta(y * b) / b, w, wake_angle) 
      + phase_(x, y + a * delta(y * b) / b, w, wake_angle) 
    ) / 2
  }

// x, y: Position in wake's coordinates
#local wake = function(x, y, w, wake_angle) {
    // Only waves for positive x, and dampen as they are further from the
    // ship's aft
    damp_half(x, 1, 0.1)
    
    // The actual waves, controlled by the (x,y) dependent phase
    * waves2_( max( 0, phase2_(x, y, w, wake_angle, 0.2 , 4) ) )

    // Dampen the waves as they get further from the main wake front
    * damp_asym(
        y
          // Less damping for larger x
          / (1 + 0.1 * x) 
              
          // Scale damping of inside wave (positive y-values) such that it 
          // always meets up in the center, even when the left and right side
          // of the wake use different settings for w and wake_angle.
          / select(y, 1, (tan(wake_angle) * x + w / cos(wake_angle))),
         
        // Damping of outside wave
        //   Aggressive damping for small x (near front of ship), and more 
        //   gradual damping further away
        1 * (1 + 50 / (1 + max(x, 0) * 1.0)), 

								// Damping of inside wave
        16
        )
  }

#local wake_mask_ = function(wx, wy, w, wake_angle) {
    select(wx, 
      1, // Plain normal in front of ship 
      select(wy, 
        // Outside wake 
        0.5 + 0.5 * min(abs(alpha_(wx, wy, w, wake_angle)) * 3, 1),
        // Inside wake
        0.5 - 0.5 * min(abs(alpha_(wx, wy, w, wake_angle)) * 9, 1)))
  }

// Transforms to the wake's coordinate system
#local x_wake = function(x, y, w, wake_angle) {
    x * cos(wake_angle) - (w - abs(y)) * sin(wake_angle)
  }

// Transforms to the wake's coordinate system
#local y_wake = function(x, y, w, wake_angle) {
    x * sin(wake_angle) + (w - abs(y)) * cos(wake_angle) 
  }

// Invokes wake function after transforming to wake coordinates.
#local wake_waves = function(x, y, w, wake_angle) {
    wake( x_wake(x, y, w, wake_angle), 
          y_wake(x, y, w, wake_angle), w, wake_angle)
  }

#local trail_angle = function(x, y, x0, w) {
    atan2(abs(y) - w, x - x0) 
  }

#local x_rotate = function(x, y, a) {
	   x * cos(a) + y * sin(a)
  }
#local y_rotate = function(x, y, a) {
	   -x * sin(a) + y * cos(a)
  }

// Creates a mask for the foam trail. 1 -> trail, 0 -> no trial
//
// x, y: Basic coordinates
// x0: Start of trail
// w: Width of half of the trail
// dy: Shift of trail, as it may not be centered around y=0
// min_angle: Angle within which trail is still at full intensity
// max_angle: Angle where trail has disappeared entirely
// angle_offset: Angle that trail makes with y-axis (it does not have to be
//               aligned with it)
#local trail_mask = 
  function(x, y, x0, w, dy, min_angle, max_angle, angle_offset) {
	   damp_sym(x - x0, 0.3) * 
    select(x - x0, 
      0,
      min(1, max(0, (trail_angle(x_rotate(x+x0, y+dy, angle_offset),
                                 y_rotate(x+x0, y+dy, angle_offset), x0, w) 
                     - max_angle) /
                    (min_angle - max_angle)
                )))
  }

#local wake_mask = function(x, y, w, wake_angle) {
    wake_mask_( x_wake(x, y, w, wake_angle), 
                y_wake(x, y, w, wake_angle), w, wake_angle)
  }

#local all_mask = function(x, y, w, wake_angle) {
	 0.2 - 0.2 * trail_mask(x, y, 1.7, 0.85, 0.10, 
	                        radians(2), radians(6), -radians(1.5))
	     + 0.8 * wake_mask(x, y, w, wake_angle) 
}

#local wind_waves2 = function(rx, ry) {
    waves2_( (rx + 0.2 * f_noise3d(rx, ry * 0.3, 0)) * 10 )
}  

#local wind_waves = function(x, y) {
    (0.2 + 0.8 * f_noise3d(y, x, 0)) *
      wind_waves2( x_rotate(x, y, radians(-20)),
                   y_rotate(x, y, radians(-20)) )
  }

#local all_waves_ = function(x, y, w, wake_angle) {
    // The wake waves
    wake_waves(x, y, w, wake_angle) +
        
    // The wind waves
    0.5 * wind_waves(x, y) * 
      // Use wake_mask to only have wind waves outside wake
      //   wake_mask < 0.55: -> 0
      //   wake_mask > 0.75: -> 1
      min(1, max( 0, (wake_mask(x, y, w, wake_angle) - 0.55) * 5) )
  }

#local all_waves = function(x, y, w, wake_angle) {
	   // Scale from [-1, 1] to [0, 1]
    0.5 + 0.5 * all_waves_(x, y, w, wake_angle)
  }

#local shadow_mask = function(x, y) {
    min(1, max(0, min(0.8 + 0.5 * x - y, 2 - y - x) ))
  }

#local all_waves_in_shadow = function(x, y, w, wake_angle) {
	   // Scale from [-1, 1] to [0, 1]
    0.5 + 0.5 * all_waves_(x, y, w, wake_angle) * shadow_mask(x, y)
  }

#local coords_ = function(x, y) {
  checker_ ( x, y, 1 )
}

#local coords_wake = function(x, y, w, wake_angle) {
  checker_ ( x_wake(x, y, w, wake_angle), 
             y_wake(x, y, w, wake_angle), 1 )
}

#local coords_innerwave = function(x, y, w, wake_angle) {
  checker_ (
    alpha_(
      x_wake(x, y, w, wake_angle),
      y_wake(x, y, w, wake_angle),
      w, wake_angle) * 18 / pi, 
    radius_(
      x_wake(x, y, w, wake_angle), 
      y_wake(x, y, w, wake_angle), 
      w, wake_angle), 1 )
}

#local ShipTransform = 
  transform {
  	 // Move center of image to origin of coordinates in main scene file.
  	 translate -x * 0.5
    translate -y * 0.25
    
    // Move center a bit to more effectively cover the scene file.
    translate -x * 2
    translate y * 1
    
    rotate z * 90
    
  }

#local f_pigment = function (x, y, w, wake_angle) {
  #if (clock < 0.33)
    all_waves(x, y, w, wake_angle)
  #else
  #if (clock < 0.67)
    all_waves_in_shadow(x, y, w, wake_angle)  
  #else
    all_mask(x, y, w, wake_angle)
  #end
  #end
  }

difference {
  plane {
    -z, 1
  }
  plane {
  	  y, 0
  }

  texture {
    pigment {
      function {
        f_pigment(x, y, 0.75, radians(19.5))
      }
    }
    finish {
      ambient 1
    }
  }
  
  transform { ShipTransform }
}


difference {
  plane {
    -z, 1
  }
  plane {
  	  -y, 0
  }

  texture {
    pigment {
      function {
        f_pigment(x, y, 0.40, radians(26.5))
      }
    }
    finish {
      ambient 1
    }
  }

  transform { ShipTransform }
}

//union {
//  box {
//    < 0,  0.75, -0.9>,
//    < 0.5, -0.25, -1.1>
//  }
//
//  box {
//    < 1,  0.75, -0.9>,
//    < 2, -1.25, -1.1>
//  }
//  
//  texture {
//    pigment { color Red }
//    finish { ambient 1 }
//  }
//
//  transform { ShipTransform }
//}


// Add a soft-border so that the borders of the image all have an RBG value
// of 50%, enabling perfect fit inside a larger entirely flat plane.
#local border = function(x, y, a, b) {
	   min(1, max(0, (Chebyshev_dist(x, y, 0, 0) - a) / (b - a)))
  }

plane {
  -z, 2
  texture {
    pigment {
      function {
        border(x, y, (Zoom/2) * 0.9, Zoom/2)
      }
      color_map {
        [0.1 color rgbt <0.5, 0.5, 0.5, 1> ]
        [0.9 White * 0.5]
      }
    }
    finish { ambient 1 }
  }
}