= #JxnPortable\docs\plot_examples\demos\WorldMap~Orthodrome.jxn

! Demonstrates the difference between the orthodrome and the loxodrome curve
! Shows both curves in equirectangular and Mercator projection and on the globe

! On the globe we can see that the orthodrome is the shortest connection though
! on the flat maps because of the distortion the loxodrome appears to be shorter


filename = "WorldMap.dig.txt"
fn = $this.getPath filename

#if @File( fn ).exists()  ! check data file

   ! read costline data:
   src = @JxnRealArrayTextFileDataSource( fn, 0, ",", " " )
   format src.readLines 5
   lon = src[1]
   lat = src[2]

#else
   javax.swing.JOptionPane.showMessageDialog( @Frame(), "Looks better with file " + filename )

   ! create geographic coordinate system (if costline is not available):
   xy = @JxnUVGrid( -90., 90., 10., -180., 180., 10., 2 )
   lat = xy.u()
   lon = xy.v()

#endif



! define slider variables:
slider = @JxnSliderPanel( $this )
latA = slider.add( "orig-lat",  51.3,  -90,  90 )
lonA = slider.add( "orig-lon",   0.0, -180, 180 )
latE = slider.add( "dest-lat",  35.5,  -90,  90 )
lonE = slider.add( "dest-lon", 139.5, -180, 180 )

! #JXN:include/P3dInclude
p3d_dz    = slider.add( "dz", 1000.,   7.1, 1000. )
p3d_alpha = slider.add( "alpha", 0.,  -90.,   90. )
p3d_phi   = slider.add( "phi",   0., -270.,  270. )
p3d = @KmgPerspectiveProjection( p3d_phi - 90, p3d_alpha, p3d_dz )



! earth radius:
re = 6.37

! Mercator projection:
mrcy = @JxnFunction( $this, "lat", "$sn = sinD lat; 0.5 log( ( 1 + $sn ) / ( 1 - $sn ) )" )
mrcp = @JxnFunction( $this, "lat", "mrcy( lat ) 180 / PI" )


! define the loxodrome:
dLon = lonE - lonA 
dLon = sw( dLon, -180, 180, dLon + 360, dLon, dLon - 360 )

dlx = max( abs( latE - latA ), abs( dLon ) )
tL = @JxnRealArrayAlgebra( 0., 1., 1/dlx )

yA = mrcy( latA )
yE = mrcy( latE )
yL = yA + ( yE - yA ) tL
latL = 2 atanD exp yL - 90

lonL = lonA + dLon tL
lonLc = clip( sw( lonL, -180, 180, lonL + 360, lonL, lonL - 360 ), -178, 178 )


! define and map a great circle (orthodrome):
tO = @JxnRealArrayAlgebra( 0., 360., 1. )
xc = re cosD tO
yc = re sinD tO

v0 = @JxnVectorAlgebra( latA, lonA, true )
v1 = @JxnVectorAlgebra( latE, lonE, true )
vn = v0 v1  ! v0 x v1
tv = @JxnVectorTransformation().rotD( 0., 90 - vn.latD(), vn.lonD()  )

xyzO = tv.map( { xc, yc, 0 tO } )

latO = asinD( xyzO[2] / re )
lonO = atan2D( xyzO[1], xyzO[0] )
lonO = clip( lonO, -179, 179 )



! draw equirectangular projection:
pf1a = plot( lon, lat, B, 0, 1. ).add( lonLc, latL, R ).add( lonO, latO, G )
pf1a.setPlotFrameTitle( "Equirectangular with orthodrome(green) and loxodrome(red)" )


! draw Mercator projection:
pf1b = plot( lon, mrcp lat, B, 0, 1. ).add( lonLc, mrcp latL, R ).add( lonO, mrcp latO, G )
pf1b.setPlotFrameTitle( "Mercator with orthodrome(green) and loxodrome(red)" )



! define mapping ( lat, lon ) -> ( x, y, z ):
xp = @JxnFunction( $this, "lon, lat", "cosD lat * cosD lon" )
yp = @JxnFunction( $this, "lon, lat", "cosD lat * sinD lon" )
zp = @JxnFunction( $this, "lon, lat", "sinD lat" )

! draw costlines:
x = re xp( lon, lat )
y = re yp( lon, lat )
z = re zp( lon, lat )
! zmin = 0.
zmin = re^2 / p3d_dz   ! <= re^2 + dh^2 = dz^2; rh^2 + zmin^2 = re^2; (dz-zmin)^2 + rh^2 = re^2
zv = sw( getZ( p3d, { x, y, z } ) - zmin, Double.NaN + z, z )  ! visible side only
pf3 = plot( p3d, 0. ).add( x, y, zv, B )
pf3.setPlotFrameTitle( "Globe with orthodrome(green) and loxodrome(red)" )

! add horizon:
! rhz = 1.
rhz = 1 / sqrt( 1 - ( re / p3d_dz )^2 )
pf3.add( curve( rhz xc, rhz yc, @Color( 204, 204, 204 ) ) )

! add orthodrome:
zOv = sw( getZ( p3d, xyzO ) - zmin, Double.NaN + tO, xyzO[2] )  ! visible -> solid line
zOh = sw( getZ( p3d, xyzO ) - zmin, xyzO[2], Double.NaN + tO )  ! hidden -> broken line
pf3.add( xyzO[0], xyzO[1], zOv, G ).add( xyzO[0], xyzO[1], zOh, G, -2 )

! add loxodrome:
xL = re xp( lonL, latL )
yL = re yp( lonL, latL )
zL = re zp( lonL, latL )
xyzL = { xL, yL, zL }
zLv = sw( getZ( p3d, xyzL ) - zmin, Double.NaN + tL, zL )  ! visible -> solid line
zLh = sw( getZ( p3d, xyzL ) - zmin, zL, Double.NaN + tL )  ! hidden -> broken line 
pf3.add( xL, yL, zLv ).add( xL, yL, zLh, R, -2 )