#!/usr/local/bin/gnuplot -persist
# set terminal canvas  rounded size 600,400 enhanced fsize 10 lw 1.6 fontscale 1 name "map_projection_3" jsdir "."
# set output 'map_projection.3.js'
unset border
set angles degrees
set key at screen 0.98, 0.95 right top vertical Right noreverse enhanced autotitle nobox
set key noinvert samplen 0.3 spacing 1 width 0 height 0 
set size ratio 0.5 1,1
set style data lines
unset xtics
unset ytics
set title "{/:Bold Albers} equal-area conic projection" 
set xrange [ * : * ] noreverse writeback
set x2range [ * : * ] noreverse writeback
set yrange [ * : * ] noreverse writeback
set y2range [ * : * ] noreverse writeback
set zrange [ * : * ] noreverse writeback
set cbrange [ * : * ] noreverse writeback
set rrange [ * : * ] noreverse writeback
set lmargin  1
set rmargin  1
set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front  noinvert bdefault
α(λ,φ) = acos(cos(φ) * cos(λ/2.))
x_W3(λ, φ) = 0.5 * (λ*cos(φ1) + (2*cos(φ)*sin(λ/2.))/sinc(α(λ,φ)))
y_W3(λ, φ) = 0.5 * (φ + sin(φ)/sinc(α(λ,φ)))
x_Hammer(λ, φ) = (2*sqrt(2)*cos(φ)*sin(λ/2)) / sqrt(1.0 + cos(φ)*cos(λ/2))
y_Hammer(λ, φ) = sqrt(2)*sin(φ) / sqrt(1.0 + cos(φ)*cos(λ/2))
Z(x,y) = sqrt( 1.0 - (x/4)**2 - (y/2)**2 )
lon(x,y) = atan( x/(ρ0-y) ) / n
lat(x,y) = asin( (C - n*n*(x**2 + (ρ0-y)**2)) / (2*n) )
θ(λ) = n*(λ-λ0)
ρ(φ) = sqrt(C - 2*n*sin(φ)) / n
x_Albers(λ, φ) = ρ(φ) * sin(θ(λ))
y_Albers(λ, φ) = ρ0 - ρ(φ)*cos(θ(λ))
clip(south) = south < -60. ? NaN : south
mouse_readout(x,y) = sprintf("Longitude %.2f Latitude %.2f", lon(x,y), lat(x,y))
π = 3.14159265358979
φ1 = 0.0
φ0 = 0.0
φ2 = 60.0
λ0 = 0.0
n = 0.433012701892219
C = 1.0
ρ0 = 2.3094010767585
## Last datafile plotted: "world.dat"
plot for [λ=-180:180:10] [φ=-60:90] '+' using (x_Albers(λ,φ)):(y_Albers(λ,φ)) with lines lc "cyan" lw .5 notitle,      for [φ = -60:90:10] [λ=-180:180:10] '+' using (x_Albers(λ,φ)):(y_Albers(λ,φ)) with lines lc "cyan" lw .5 notitle,      'world.dat' using (x_Albers($1,clip($2))):(y_Albers($1,clip($2))) with filledcurve fc "brown" title "fill",      'world.dat' using (x_Albers($1,clip($2))):(y_Albers($1,clip($2))) with lines lc "black" title "outline"