Generated: Sat Oct 24 16:35:12 2020 from cmpdist.pl 2020/05/18 4.8 KB. text copy
#!/usr/bin/perl -w # NAME: distance.pl # AIM: Use perl trib function for distance London to Tokyo... # from : http://perldoc.perl.org/Math/Trig.html use strict; use warnings; use Math::Trig qw(great_circle_distance great_circle_direction deg2rad rad2deg); my $os = $^O; my $perl_dir = '/home/geoff/bin'; my $PATH_SEP = '/'; my $temp_dir = '/tmp'; if ($os =~ /win/i) { $perl_dir = 'C:\GTools\perl'; $temp_dir = $perl_dir; $PATH_SEP = "\\"; } unshift(@INC, $perl_dir); require 'fg_wsg84.pl' or die "Unable to load fg_wsg84.pl ...\n"; # 1 knots = 1.85200 kph my $K2KPH = 1.85200; my $lonlon = -0.5; my $lonlat = 51.3; my $toklon = 139.8; my $toklat = 35.7; my $g_ias = 80; # Knots - c152-c172 my $g_speed = $g_ias * $K2KPH; # Kilometers/Hour my $CP_EPSILON = 0.00001; # Notice the 90 - latitude: phi zero is at the North Pole sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) } sub show_sg_distance_vs_est($$$$$$) { my ($lat1,$lat2,$lon1,$lon2,$estdist,$esthdg) = @_; my ($sg_az1,$sg_az2,$sg_dist,$res,$sg_km,$sg_ikm); my ($eddiff,$edpc,$ahdg,$ahdiff,$hdgpc,$adhdg,$distmsg,$hdgmsg); $res = fg_geo_inverse_wgs_84 ($lat1,$lon1,$lat2,$lon2,\$sg_az1,\$sg_az2,\$sg_dist); $sg_km = $sg_dist / 1000; $sg_ikm = int($sg_km + 0.5); $sg_az1 = int(($sg_az1 * 10) + 0.05) / 10; $sg_az2 = int(($sg_az2 * 10) + 0.05) / 10; $eddiff = abs($sg_km - $estdist); $edpc = ($eddiff / $sg_km) * 100; $edpc = int(($edpc * 10) + 0.05) / 10; if ($eddiff < $CP_EPSILON) { $edpc = "SAME"; } elsif ($edpc < 1) { $edpc = 'lt 1'; } $edpc = " $edpc" while (length($edpc) < 5); $adhdg = abs($sg_az1 - $esthdg); $ahdiff = ($adhdg / $sg_az1); $hdgpc = ($ahdiff / 360); $hdgpc = int(($hdgpc * 10) + 0.05) / 10; if ($ahdiff < $CP_EPSILON) { $hdgpc = "SAME"; } elsif ($hdgpc < 1) { $hdgpc = 'lt 1'; } $hdgpc = " $hdgpc" while (length($hdgpc) < 5); print "SG_Dist : $sg_ikm kilometers ($sg_km) ($edpc \% diff)\n"; print "SG_Head : $sg_az1 degs (inverse $sg_az2) ($hdgpc \% diff)\n"; } sub show_distance($$$$) { my ($lon1,$lat1,$lon2,$lat2) = @_; my @Pos1 = NESW( $lon1, $lat1); my @Pos2 = NESW( $lon2, $lat2); my $d_km = great_circle_distance(@Pos1, @Pos2, 6378); # About 9600 km. my $t_degs = rad2deg(great_circle_distance(@Pos1, @Pos2)); my $hdg = rad2deg(great_circle_direction(@Pos1, @Pos2)); $hdg = int(($hdg * 10) + 0.05) / 10; my $hrs = $d_km / $g_speed; my $eta = get_hhmmss($hrs); my $ikm = int($d_km + 0.5); #$degs = int($degs + 0.5); #$degs = int($degs + 0.5); $t_degs = sprintf("%0.6f",$t_degs); print "From (lon,lat) $lon1,$lon2 to $lon2,$lat2 is about -\n"; print "Distance: $ikm kilometers ($d_km)\n"; print "Heading : $hdg, for $t_degs degs\n"; print "ETA : $eta, at $g_ias Knots\n"; # print "km $km / spd $g_speed = $hours\n"; show_sg_distance_vs_est($lat1,$lat2,$lon1,$lon2,$d_km,$t_degs); # ================================================= } sub get_hhmmss($) { my ($hours) = @_; my $hrs = int($hours); my $mins = ($hours * 60) - ($hrs * 60); my $min = int($mins); my $secs = ($mins * 60) - ($min * 60); my $days = 0; if ($hrs >= 24) { $days++; $hrs -= 24; } $min = ($min < 10) ? "0$min" : $min; $secs = int( $secs + 0.5 ); $secs = ($secs < 10) ? "0$secs" : $secs; $hrs = ($hrs < 10) ? "0$hrs" : $hrs; my $stg = ''; $stg .= "Days $days, " if ($days); $stg .= "$hrs:$min:$secs"; return $stg; } my ($rkm,$heading,@L,@T,$km,$degs,$argc,$hours,$hms); $argc = scalar @ARGV; if ($argc == 4) { $lonlon = $ARGV[0]; $lonlat = $ARGV[1]; $toklon = $ARGV[2]; $toklat = $ARGV[3]; @L = NESW( $lonlon, $lonlat); @T = NESW( $toklon, $toklat); $km = great_circle_distance(@L, @T, 6378); # About 9600 km. $degs = rad2deg(great_circle_distance(@L, @T)); $heading = rad2deg(great_circle_direction(@L, @T)); $heading = int(($heading * 10) + 0.05) / 10; $hours = $km / $g_speed; $hms = get_hhmmss($hours); $rkm = int($km + 0.5); #$degs = int($degs + 0.5); #$degs = int($degs + 0.5); $degs = sprintf("%0.6f",$degs); print "From (lon,lat) $lonlon,$lonlat to $toklon,$toklat is about -\n"; print "Distance: $rkm kilometers ($km)\n"; print "Heading : $heading, for $degs degs\n"; print "ETA : $hms, at $g_ias Knots\n"; # print "km $km / spd $g_speed = $hours\n"; } else { print "ERROR: Need 4 arguments, not $argc!\n" if ($argc); print "Given two pairs of longitude, latitude, and will calculate the distance, and heading...\n"; show_distance( $lonlon, $lonlat, $toklon, $toklat ); } # eof - distance.pl