#!/usr/bin/perl -w
# NAME: magvar.pl
# AIM: Given lat, lon, hh dd mm yy output magnetic variation
# 09/10/2014 - Attempt to fix 
# Use of uninitialized value in multiplication (*) at magvar.pl line 335.
# Use of uninitialized value in multiplication (*) at magvar.pl line 338.
# 13/11/2013 geoff mclane http://geoffair.net/mperl
use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Data::Dumper;
use Cwd;
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 'lib_utils.pl' or die "Unable to load 'lib_utils.pl' Check paths in \@INC...\n";
# log file stuff
our ($LF);
my $pgmname = $0;
if ($pgmname =~ /(\\|\/)/) {
    my @tmpsp = split(/(\\|\/)/,$pgmname);
    $pgmname = $tmpsp[-1];
}
my $outfile = $temp_dir.$PATH_SEP."temp.$pgmname.txt";
open_log($outfile);

# user variables
my $VERS = "0.0.1 2013-03-17";
my $load_log = 0;
my $in_file = '';
my $verbosity = 0;
my $out_file = '';

# ### DEBUG ###
my $debug_on = 0;
my $def_file = 'def_file';

### program variables
my @warnings = ();
my $cwd = cwd();

sub VERB1() { return $verbosity >= 1; }
sub VERB2() { return $verbosity >= 2; }
sub VERB5() { return $verbosity >= 5; }
sub VERB9() { return $verbosity >= 9; }

sub show_warnings($) {
    my ($val) = @_;
    if (@warnings) {
        prt( "\nGot ".scalar @warnings." WARNINGS...\n" );
        foreach my $itm (@warnings) {
           prt("$itm\n");
        }
        prt("\n");
    } else {
        prt( "\nNo warnings issued.\n\n" ) if (VERB9());
    }
}

sub pgm_exit($$) {
    my ($val,$msg) = @_;
    if (length($msg)) {
        $msg .= "\n" if (!($msg =~ /\n$/));
        prt($msg);
    }
    show_warnings($val);
    close_log($outfile,$load_log);
    exit($val);
}


sub prtw($) {
   my ($tx) = shift;
   $tx =~ s/\n$//;
   prt("$tx\n");
   push(@warnings,$tx);
}

sub process_in_file($) {
    my ($inf) = @_;
    if (! open INF, "<$inf") {
        pgm_exit(1,"ERROR: Unable to open file [$inf]\n"); 
    }
    my @lines = <INF>;
    close INF;
    my $lncnt = scalar @lines;
    prt("Processing $lncnt lines, from [$inf]...\n");
    my ($line,$inc,$lnn);
    $lnn = 0;
    foreach $line (@lines) {
        chomp $line;
        $lnn++;
        if ($line =~ /\s*#\s*include\s+(.+)$/) {
            $inc = $1;
            prt("$lnn: $inc\n");
        }
    }
}

my $lat_deg = 9999;
my $lon_deg = 9999;
my $hour = 0;
my $day2 = 13;
my $month2 = 11;
my $year2 = 13;  ### 2013;
my $day = 10;
my $month = 10;
my $year = 14;  ### 2013;
#/* Convert date to Julian day    1950-2049 */
sub yymmdd_to_julian_days($$$) {
    my ($yy, $mm, $dd) = @_;
    my $jd;
 
    $yy = ($yy < 50) ? (2000 + $yy) : (1900 + $yy);
    $jd = $dd - 32075 + 1461 * ($yy + 4800 + ($mm - 14) / 12 ) / 4;
    $jd = $jd + 367 * ($mm - 2 - ($mm - 14) / 12*12) / 12;
    $jd = $jd - 3 * (($yy + 4900 + ($mm - 14) / 12) / 100) / 4;
    #/* printf("julian date = %d\n", jd ); */
    return $jd;
} 

my $SGD_PI = 3.1415926535;
my $SGD_DEGREES_TO_RADIANS = $SGD_PI / 180.0;
my $SGD_RADIANS_TO_DEGREES = 180.0 / $SGD_PI;
#var = calc_magvar( SGD_DEGREES_TO_RADIANS * lat_deg, SGD_DEGREES_TO_RADIANS * lon_deg, h,
#		   yymmdd_to_julian_days(yy,mm,dd), field );

my $nmax = 12;
my $a = 6378.137;       #/* semi-major axis (equatorial radius) of WGS84 ellipsoid */
my $b = 6356.7523142;   #/* semi-minor axis referenced to the WGS84 ellipsoid */
my $r_0 = 6371.2;	    #/* standard Earth magnetic reference radius  */

my @gnm_wmm2005 =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-29556.8, -1671.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-2340.6, 3046.9, 1657.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [1335.4, -2305.1, 1246.7, 674.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [919.8, 798.1, 211.3, -379.4, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-227.4, 354.6, 208.7, -136.5, -168.3, -14.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [73.2, 69.7, 76.7, -151.2, -14.9, 14.6, -86.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [80.1, -74.5, -1.4, 38.5, 12.4, 9.5, 5.7, 1.8, 0.0, 0.0, 0.0, 0.0, 0.0],
    [24.9, 7.7, -11.6, -6.9, -18.2, 10.0, 9.2, -11.6, -5.2, 0.0, 0.0, 0.0, 0.0],
    [5.6, 9.9, 3.5, -7.0, 5.1, -10.8, -1.3, 8.8, -6.7, -9.1, 0.0, 0.0, 0.0],
    [-2.3, -6.3, 1.6, -2.6, 0.0, 3.1, 0.4, 2.1, 3.9, -0.1, -2.3, 0.0, 0.0],
    [2.8, -1.6, -1.7, 1.7, -0.1, 0.1, -0.7, 0.7, 1.8, 0.0, 1.1, 4.1, 0.0],
    [-2.4, -0.4, 0.2, 0.8, -0.3, 1.1, -0.5, 0.4, -0.3, -0.3, -0.1, -0.3, -0.1]
);

my @hnm_wmm2005 =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 5079.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -2594.7, -516.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -199.9, 269.3, -524.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 281.5, -226.0, 145.8, -304.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 42.4, 179.8, -123.0, -19.5, 103.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -20.3, 54.7, 63.6, -63.4, -0.1, 50.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -61.5, -22.4, 7.2, 25.4, 11.0, -26.4, -5.1, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 11.2, -21.0, 9.6, -19.8, 16.1, 7.7, -12.9, -0.2, 0.0, 0.0, 0.0, 0.0],
    [0.0, -20.1, 12.9, 12.6, -6.7, -8.1, 8.0, 2.9, -7.9, 6.0, 0.0, 0.0, 0.0],
    [0.0, 2.4, 0.2, 4.4, 4.8, -6.5, -1.1, -3.4, -0.8, -2.3, -7.9, 0.0, 0.0],
    [0.0, 0.3, 1.2, -0.8, -2.5, 0.9, -0.6, -2.7, -0.9, -1.3, -2.0, -1.2, 0.0],
    [0.0, -0.4, 0.3, 2.4, -2.6, 0.6, 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8]
);

my @gtnm_wmm2005 =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [8.0, 10.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-15.1, -7.8, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.4, -2.6, -1.2, -6.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-2.5, 2.8, -7.0, 6.2, -3.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-2.8, 0.7, -3.2, -1.1, 0.1, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-0.7, 0.4, -0.3, 2.3, -2.1, -0.6, 1.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.2, -0.1, -0.3, 1.1, 0.6, 0.5, -0.4, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.1, 0.3, -0.4, 0.3, -0.3, 0.2, 0.4, -0.7, 0.4, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
);

my @htnm_wmm2005 =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -20.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -23.2, -14.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 5.0, -7.0, -0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 2.2, 1.6, 5.8, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 1.7, 2.1, 4.8, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -0.6, -1.9, -0.4, -0.5, -0.3, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.6, 0.4, 0.2, 0.3, -0.8, -0.2, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, -0.2, 0.1, 0.3, 0.4, 0.1, -0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
);



my @P =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
);
my @DP =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
);
my @gnm =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
);

my @hnm =
(
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
);

my @sm = ([1,2,3,4,5,6,7,8,9,10,11,12,13]);
my @cm = ([1,2,3,4,5,6,7,8,9,10,11,12,13]);
my @root = ([1,2,3,4,5,6,7,8,9,10,11,12,13]);
my @roots = ();

sub create_roots_array() {
    my ($i,$j,$k,$n);
    prt("Create roots array for nmax $nmax\n");
    for ($i = 0; $i <= $nmax; $i++) {
        for ($j = 0; $j <= $nmax; $j++) {
            for ($k = 0; $k < 2; $k++) {
                $roots[$i][$j][$k] = 1;
            }
        }
    }
    for ($i = 0; $i <= $nmax; $i++) {
        my $ii = $i*$i;
        ##$j = ((($i + 1) > 2) ? ($i + 1) : 2);
        for ($j = 0; $j <= $nmax; $j++) {
            $n = (($i + $j + 1) > 2) ? ($i + $j + 1) : 2;
            my $s1 = sqrt(($n-1)*($n-1) - $ii);
            my $s2 = (1.0 / sqrt( $n*$n - $ii));
            $roots[$i][$j][0] = $s1;
            $roots[$i][$j][1] = $s2;
        }
    }
}

#/*
# * return variation (in radians) given geodetic latitude (radians),
# * longitude(radians), height (km) and (Julian) date
# * N and E lat and long are positive, S and W negative
#*/

sub calc_magvar($$$$) {
    my ($lat, $lon, $h, $dat ) = @_;
    #/* output field B_r,B_th,B_phi,B_x,B_y,B_z */
    my ($n,$m);
    #/* reference date for current model is 1 januari 2005 */
    my $date0_wmm2005 = yymmdd_to_julian_days(5,1,1);

    my ($yearfrac,$sr,$r,$theta,$c,$s,$psi,$fn,$fn_0,$B_r,$B_theta,$B_phi,$X,$Y,$Z);
    my ($sinpsi, $cospsi, $inv_s);
    my ($msg);
    #static int been_here = 0;

    my $sinlat = sin($lat);
    my $coslat = cos($lat);

    #/* convert to geocentric coords: */
    #// sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
    $sr = sqrt($a*$a*$coslat*$coslat + $b*$b*$sinlat*$sinlat);
    #/* sr is effective radius */
    $theta = atan2($coslat * ($h*$sr + $a*$a),
		  $sinlat * ($h*$sr + $b*$b));
    #/* theta is geocentric co-latitude */

    $r = $h*$h + 2.0*$h * $sr +
	($a*$a*$a*$a - ( $a*$a*$a*$a - $b*$b*$b*$b ) * $sinlat*$sinlat ) / 
	($a*$a - ($a*$a - $b*$b) * $sinlat*$sinlat );

    $r = sqrt($r);

    #/* r is geocentric radial distance */
    $c = cos($theta);
    $s = sin($theta);
    #/* protect against zero divide at geographic poles */
    $inv_s =  1.0 / ($s + ($s == 0.0)*1.0e-8); 

    #/* zero out arrays */
    for ( $n = 0; $n <= $nmax; $n++ ) {
    	### for ( $m = 0; $m <= $n; $m++ ) {  # CHECK ME
    	for ( $m = 0; $m <= $nmax; $m++ ) {
	        $P[$n][$m] = 0;
	        $DP[$n][$m] = 0;
            #if (($m == 0) && ($n == 2)) {
            #    prt("For m=$m, n=$n set DP[n][m] = 0\n");
            #}
	    }
    }

    #/* diagonal elements */
    $P[0][0] = 1;
    $P[1][0] = $c;
    $P[1][1] = $s;

    $DP[0][0] = 0;
    $DP[1][0] = -$s;
    $DP[1][1] = $c;

    #// these values will not change for subsequent function calls
    #if( !been_here ) {
        for ( $n = 2; $n <= $nmax; $n++ ) {
            $root[$n] = sqrt((2.0*$n-1) / (2.0*$n));
        }

#        for ( $m = 0; $m <= $nmax; $m++ ) {
#            my $mm = $m*$m;
#            my ($s1,$s2);
#            $n = ((($m + 1) > 2) ? ($m + 1) : 2);
#            ##for ( $n = SG_MAX2($m + 1, 2); $n <= $nmax; $n++ ) {
#            for ( ; $n <= $nmax; $n++ ) {
#                $s1 = sqrt(($n-1)*($n-1) - $mm);
#                $s2 = (1.0 / sqrt( $n*$n - $mm));
#                $roots[$m][$n][0] = $s1;
#                $roots[$m][$n][1] = $s2;
#            }
#        }
	#been_here = 1;
    #}

    for ( $n = 2; $n <= $nmax; $n++ ) {
    	#// double root = sqrt((2.0*n-1) / (2.0*n));
	    $P[$n][$n] = $P[$n-1][$n-1] * $s * $root[$n];
	    $DP[$n][$n] = ($DP[$n-1][$n-1] * $s + $P[$n-1][$n-1] * $c) * $root[$n];
    }

    #/* lower triangle */
    for ( $m = 0; $m <= $nmax; $m++ ) {
        #// double mm = m*m;
        $n = ((($m + 1) > 2) ? $m + 1 : 2);
        #for ( $n = SG_MAX2($m + 1, 2); $n <= $nmax; n++ ) {
        for ( ; $n <= $nmax; $n++ ) {
            #// double root1 = sqrt((n-1)*(n-1) - mm);
            #// double root2 = 1.0 / sqrt( n*n - mm);
            my $rv0 = $roots[$m][$n][0];
            my $rv1 = $roots[$m][$n][1];
            my $p1 = $P[$n-1][$m];
            my $p2 = $P[$n-2][$m];
            my $dp1 = $DP[$n-1][$m];
            my $dp2 = $DP[$n-2][$m];
            if (defined $rv0 && defined $rv1 && defined $p1 && defined $p2 && defined $dp1 && defined $dp2) {
                ### prt("m=$m, n=$n, rv0=$rv0, rv1=$rv1, c=$c, s=$s\n");
                $P[$n][$m] = ($P[$n-1][$m] * $c * (2.0*$n-1) -
                       $P[$n-2][$m] * $rv0) * $rv1;

                $DP[$n][$m] = (($DP[$n-1][$m] * $c - $P[$n-1][$m] * $s) *
                    (2.0*$n-1) - $DP[$n-2][$m] * $rv0) * $rv1;
            } else {
                $msg = '';
                $msg .= "no rv0 " if ( !(defined $rv0) );
                $msg .= "no rv1 " if ( !(defined $rv1) );
                $msg .= "no p1 " if ( !(defined $p1) );
                $msg .= "no p2 " if ( !(defined $p2) );
                $msg .= "no dp1 " if ( !(defined $dp1) );
                $msg .= "no dp2 " if ( !(defined $dp2) );
                pgm_exit(1,"Error: NOT set for m=$m, n=$n $msg\n");
            }
        }
    }

    #/* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
    #   achieved by adjusting the coefficients at time t0 for linear secular variation */
    #/* WMM2005 */
    $yearfrac = ($dat - $date0_wmm2005) / 365.25;
    for ( $n = 1; $n <= $nmax; $n++ ) {
	for ( $m = 0; $m <= $nmax; $m++ ) {
	    $gnm[$n][$m] = $gnm_wmm2005[$n][$m] + $yearfrac * $gtnm_wmm2005[$n][$m];
	    $hnm[$n][$m] = $hnm_wmm2005[$n][$m] + $yearfrac * $htnm_wmm2005[$n][$m];
	}
    }

    #/* compute sm (sin(m lon) and cm (cos(m lon)) */
    for ( $m = 0; $m <= $nmax; $m++ ) {
    	$sm[$m] = sin($m * $lon);
	    $cm[$m] = cos($m * $lon);
    }

    #/* compute B fields */
    $B_r = 0.0;
    $B_theta = 0.0;
    $B_phi = 0.0;
    $fn_0 = $r_0/$r;
    $fn = $fn_0 * $fn_0;

    for ( $n = 1; $n <= $nmax; $n++ ) {
        my $c1_n=0;
        my $c2_n=0;
        my $c3_n=0;
        for ( $m = 0; $m <= $n; $m++ ) {
            my $tmp = ($gnm[$n][$m] * $cm[$m] + $hnm[$n][$m] * $sm[$m]); 
            $c1_n=$c1_n + $tmp * $P[$n][$m];
            $c2_n=$c2_n + $tmp * $DP[$n][$m];
            $c3_n=$c3_n + $m * ($gnm[$n][$m] * $sm[$m] - $hnm[$n][$m] * $cm[$m]) * $P[$n][$m];
        }
        #// fn=pow(r_0/r,n+2.0);
        $fn *= $fn_0;
        $B_r = $B_r + ($n + 1) * $c1_n * $fn;
        $B_theta = $B_theta - $c2_n * $fn;
        $B_phi = $B_phi + $c3_n * $fn * $inv_s;
    }

    #/* Find geodetic field components: */
    $psi = $theta - (($SGD_PI / 2.0) - $lat);
    $sinpsi = sin($psi);
    $cospsi = cos($psi);
    $X = -$B_theta * $cospsi - $B_r * $sinpsi;
    $Y = $B_phi;
    $Z = $B_theta * $sinpsi - $B_r * $cospsi;

    #field[0]=B_r;
    #field[1]=B_theta;
    #field[2]=B_phi;
    #field[3]=X;
    #field[4]=Y;
    #field[5]=Z;   /* output fields */

    #/* find variation in radians */
    #/* return zero variation at magnetic pole X=Y=0. */
    #/* E is positive */
    return ($X != 0. || $Y != 0.0) ? atan2($Y, $X) : 0.0;  
}

########################################################################################
### Another try to convert C++ code to perl
#/*
# * return variation (in radians) given geodetic latitude (radians),
# * longitude(radians), height (km) and (Julian) date
# * N and E lat and long are positive, S and W negative
#*/
# double calc_magvar( double lat, double lon, double h, long dat, double* field )
my $been_here = 0;

sub SG_MAX2($$) {
    my ($a,$b) = @_;
    return ($a > $b) ? $a : $b;
}

sub calc_magvar2($$$$$) {
    my ($lat, $lon, $h, $dat, $field) = @_;
    #/* output field B_r,B_th,B_phi,B_x,B_y,B_z */
    my ($n,$m);
    #/* reference date for current model is 1 januari 2005 */
    my $date0_wmm2005 = yymmdd_to_julian_days(5,1,1);

    my ($yearfrac,$sr,$r,$theta,$c,$s,$psi,$fn,$fn_0,$B_r,$B_theta,$B_phi,$X,$Y,$Z);
    my ($sinpsi, $cospsi, $inv_s);


    my $sinlat = sin($lat);
    my $coslat = cos($lat);

    #/* convert to geocentric coords: */
    #// sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
    $sr = sqrt($a*$a*$coslat*$coslat + $b*$b*$sinlat*$sinlat);
    #/* sr is effective radius */
    $theta = atan2($coslat * ($h*$sr + $a*$a),
		  $sinlat * ($h*$sr + $b*$b));
    #/* theta is geocentric co-latitude */

    $r = $h*$h + 2.0*$h * $sr +
	    ($a*$a*$a*$a - ( $a*$a*$a*$a - $b*$b*$b*$b ) * $sinlat*$sinlat ) / 
	    ($a*$a - ($a*$a - $b*$b) * $sinlat*$sinlat );

    $r = sqrt($r);

    #/* r is geocentric radial distance */
    $c = cos($theta);
    $s = sin($theta);
    #/* protect against zero divide at geographic poles */
    $inv_s =  1.0 / ($s + ($s == 0.)*1.0e-8); 

    #/* zero out arrays */
    for ( $n = 0; $n <= $nmax; $n++ ) {
        for ( $m = 0; $m <= $n; $m++ ) {
            $P[$n][$m] = 0;
            $DP[$n][$m] = 0;
        }
    }

    #/* diagonal elements */
    $P[0][0] = 1;
    $P[1][1] = $s;
    $DP[0][0] = 0;
    $DP[1][1] = $c;
    $P[1][0] = $c ;
    $DP[1][0] = -$s;

    #// these values will not change for subsequent function calls
    if( !$been_here ) {
        for ( $n = 2; $n <= $nmax; $n++ ) {
            $root[$n] = sqrt((2.0*$n-1) / (2.0*$n));
        }

        for ( $m = 0; $m <= $nmax; $m++ ) {
            my $mm = $m*$m;
            for ( $n = SG_MAX2($m + 1, 2); $n <= $nmax; $n++ ) {
                $roots[$m][$n][0] = sqrt(($n-1)*($n-1) - $mm);
                $roots[$m][$n][1] = 1.0 / sqrt( $n*$n - $mm);
            }
        }
        $been_here = 1;
    }

    for ( $n=2; $n <= $nmax; $n++ ) {
    	#// double root = sqrt((2.0*n-1) / (2.0*n));
	    $P[$n][$n] = $P[$n-1][$n-1] * $s * $root[$n];
	    $DP[$n][$n] = ($DP[$n-1][$n-1] * $s + $P[$n-1][$n-1] * $c) *
	        $root[$n];
    }

    #/* lower triangle */
    for ( $m = 0; $m <= $nmax; $m++ ) {
        #// double mm = m*m;
        for ( $n = SG_MAX2($m + 1, 2); $n <= $nmax; $n++ ) {
            #// double root1 = sqrt((n-1)*(n-1) - mm);
            #// double root2 = 1.0 / sqrt( n*n - mm);
            $P[$n][$m] = ($P[$n-1][$m] * $c * (2.0*$n-1) -
                   $P[$n-2][$m] * $roots[$m][$n][0]) *
            $roots[$m][$n][1];

            $DP[$n][$m] = (($DP[$n-1][$m] * $c - $P[$n-1][$m] * $s) *
                (2.0*$n-1) - $DP[$n-2][$m] * $roots[$m][$n][0]) *
            $roots[$m][$n][1];
        }
    }

    #/* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
    #   achieved by adjusting the coefficients at time t0 for linear secular variation */
    #/* WMM2005 */
    $yearfrac = ($dat - $date0_wmm2005) / 365.25;
    for ( $n = 1; $n <= $nmax; $n++ ) {
        for ( $m = 0; $m <= $nmax; $m++ ) {
            $gnm[$n][$m] = $gnm_wmm2005[$n][$m] + $yearfrac * $gtnm_wmm2005[$n][$m];
            $hnm[$n][$m] = $hnm_wmm2005[$n][$m] + $yearfrac * $htnm_wmm2005[$n][$m];
        }
    }

    #/* compute sm (sin(m lon) and cm (cos(m lon)) */
    for ( $m = 0; $m <= $nmax; $m++ ) {
        $sm[$m] = sin($m * $lon);
        $cm[$m] = cos($m * $lon);
    }

    #/* compute B fields */
    $B_r = 0.0;
    $B_theta = 0.0;
    $B_phi = 0.0;
    $fn_0 = $r_0/$r;
    $fn = $fn_0 * $fn_0;

    for ( $n = 1; $n <= $nmax; $n++ ) {
        my $c1_n=0;
        my $c2_n=0;
        my $c3_n=0;
        for ( $m = 0; $m <= $n; $m++ ) {
            my $tmp = ($gnm[$n][$m] * $cm[$m] + $hnm[$n][$m] * $sm[$m]); 
            $c1_n=$c1_n + $tmp * $P[$n][$m];
            $c2_n=$c2_n + $tmp * $DP[$n][$m];
            $c3_n=$c3_n + $m * ($gnm[$n][$m] * $sm[$m] - $hnm[$n][$m] * $cm[$m]) * $P[$n][$m];
        }
        #// fn=pow(r_0/r,n+2.0);
        $fn *= $fn_0;
        $B_r = $B_r + ($n + 1) * $c1_n * $fn;
        $B_theta = $B_theta - $c2_n * $fn;
        $B_phi = $B_phi + $c3_n * $fn * $inv_s;
    }

    #/* Find geodetic field components: */
    $psi = $theta - (($SGD_PI / 2.0) - $lat);
    $sinpsi = sin($psi);
    $cospsi = cos($psi);
    $X = -$B_theta * $cospsi - $B_r * $sinpsi;
    $Y = $B_phi;
    $Z = $B_theta * $sinpsi - $B_r * $cospsi;

    ${$field}[0] = $B_r;
    ${$field}[1] = $B_theta;
    ${$field}[2] = $B_phi;
    ${$field}[3] = $X;
    ${$field}[4] = $Y;
    ${$field}[5] = $Z;   #/* output fields */

    #/* find variation in radians */
    #/* return zero variation at magnetic pole X=Y=0. */
    #/* E is positive */
    return ($X != 0 || $Y != 0) ? atan2($Y,$X) : 0;  
}


########################################################################################



sub do_test() {
    create_roots_array();
    ##prt(Dumper(\@roots));
    ##exit(0);
    my $lat = 48.7269692925;
    my $lon = 2.3699923175;
    my $alt = 0;    # km
    my $yy = ($year < 50) ? (2000 + $year) : (1900 + $year);

    my $jdate = yymmdd_to_julian_days($year, $month, $day);
    my $var = calc_magvar( $SGD_DEGREES_TO_RADIANS * 48.7269692925,
        $SGD_DEGREES_TO_RADIANS * 2.3699923175, 0, $jdate );

    prt("For lat,lon,alt $lat,$lon,$alt, date $yy/$month/$day, jdata $jdate,\nvar = ".$SGD_RADIANS_TO_DEGREES * $var." degrees\n");
}

sub do_test2() {
    ##prt(Dumper(\@roots));
    ##exit(0);
    my $lat = 48.7269692925;
    my $lon = 2.3699923175;
    my $alt = 0;    # km
    my @fields = ();
    my $yy = ($year < 50) ? (2000 + $year) : (1900 + $year);
    my $jdate = yymmdd_to_julian_days($year, $month, $day);
    my ($i,$len,$var);

    $var = calc_magvar2( $SGD_DEGREES_TO_RADIANS * $lat,
        $SGD_DEGREES_TO_RADIANS * $lon, $alt, $jdate, \@fields );
    prt("For lat,lon,alt $lat,$lon,$alt, date $yy/$month/$day, jdata $jdate,\nvar = ".$SGD_RADIANS_TO_DEGREES * $var." degrees\n");
    $len = scalar @fields;
    prt("$len fields: ");
    for ($i = 0; $i < $len; $i++) {
        prt(" ".$fields[$i]);
    }
    prt("\n");

    $lat = 37.618674211;
    $lon = -122.375007609;
    $var = calc_magvar2( $SGD_DEGREES_TO_RADIANS * $lat,
        $SGD_DEGREES_TO_RADIANS * $lon, $alt, $jdate, \@fields );
    prt("For lat,lon,alt $lat,$lon,$alt, date $yy/$month/$day, jdata $jdate,\nvar = ".$SGD_RADIANS_TO_DEGREES * $var." degrees\n");
    $len = scalar @fields;
    prt("$len fields: ");
    for ($i = 0; $i < $len; $i++) {
        prt(" ".$fields[$i]);
    }
    prt("\n");
}

#########################################
### MAIN ###
#parse_args(@ARGV);
#process_in_file($in_file);
#do_test();
do_test2();
pgm_exit(0,"");
########################################

sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have a following argument!\n") if (!@av);
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg);
    while (@av) {
        $arg = $av[0];
        if ($arg =~ /^-/) {
            $sarg = substr($arg,1);
            $sarg = substr($sarg,1) while ($sarg =~ /^-/);
            if (($sarg =~ /^h/i)||($sarg eq '?')) {
                give_help();
                pgm_exit(0,"Help exit(0)");
            } elsif ($sarg =~ /^v/) {
                if ($sarg =~ /^v.*(\d+)$/) {
                    $verbosity = $1;
                } else {
                    while ($sarg =~ /^v/) {
                        $verbosity++;
                        $sarg = substr($sarg,1);
                    }
                }
                prt("Verbosity = $verbosity\n") if (VERB1());
            } elsif ($sarg =~ /^l/) {
                if ($sarg =~ /^ll/) {
                    $load_log = 2;
                } else {
                    $load_log = 1;
                }
                prt("Set to load log at end. ($load_log)\n") if (VERB1());
            } elsif ($sarg =~ /^o/) {
                need_arg(@av);
                shift @av;
                $sarg = $av[0];
                $out_file = $sarg;
                prt("Set out file to [$out_file].\n") if (VERB1());
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } else {
            $in_file = $arg;
            prt("Set input to [$in_file]\n") if (VERB1());
        }
        shift @av;
    }

    if ($debug_on) {
        prtw("WARNING: DEBUG is ON!\n");
        if ((length($in_file) ==  0) && $debug_on) {
            $in_file = $def_file;
            prt("Set DEFAULT input to [$in_file]\n");
        }
    }
    if (length($in_file) ==  0) {
        pgm_exit(1,"ERROR: No input files found in command!\n");
    }
    if (! -f $in_file) {
        pgm_exit(1,"ERROR: Unable to find in file [$in_file]! Check name, location...\n");
    }
}

sub give_help {
    prt("$pgmname: version $VERS\n");
    prt("Usage: $pgmname [options] in-file\n");
    prt("Options:\n");
    prt(" --help  (-h or -?) = This help, and exit 0.\n");
    prt(" --verb[n]     (-v) = Bump [or set] verbosity. def=$verbosity\n");
    prt(" --load        (-l) = Load LOG at end. ($outfile)\n");
    prt(" --out <file>  (-o) = Write output to this file.\n");
}

# eof - template.pl