#!/usr/bin/perl -w
use diagnostics;
use English;
# r.downscale.pl - GRASS GIS script that sets up r.mapcalc to downscale 
# resolution for raster variables such as average temperature or degree-days,
# that tend to be highly correlated to elevation. 
#
# example use: you have low resolution (4 km) avg temperature and elev. layers,
#              and a high res. elev. (100m) layer. On the first pass, build the 
#              regression function for the low res. temp. (y) vs. elev. (x) using wt.d 
#              linear regression (5 x 5 cell local areas). On the second
#              pass, after changing to the higher resolution, apply the model
#              to the higher res elev data, using a weighted average
#              over 5 x 5 cell local areas to provide a smoothing. Wt. values are
#              selectable for both passes, allowing one to emulate 1/distance,
#              1/distance2 and other relationships.
#
# This perl program selects and sends an equation to r.mapcalc, by:
#	1. remove comments: s/#.*$// (or s/#.*?$//mg, if it's in one string)
#       2. drop blank lines; r.mapcalc treats them as EOFs
#	3. substitute input parameters: print in double-quotes
#       4. feed to r.mapcalc's stdin
#
# Developed and released to the public under a GNU public license by 
# Leonard Coop, Oregon State University, coopl@bcc.orst.edu
# Perl to r.mapcalc script orig. by Dan Upper
#
$allparams=<<END;
  ######################################################################
  ##                    CHANGE THE DEFAULTS HERE                      ##
  ##    This is not perl code -- it's key=value pairs, parsed below   ##
  ##                    Use 'NULL' in place of ""                     ##
  ######################################################################
  #
  #parameter =default    # description (all filenames are GRASS raster maplayers) 
  #--------------------------------------------------------------------------------
  calc       =BUILD      # req: function performed (BUILD  or EST)
  input      =H_50       # req: BUILD low res filename to be downscaled
  elev_lo    =usdem      # req: BUILD low res elev filename
  elev_hi    =DEM30m     # req: EST high res elev filename
  output     =HIDD5d     # req: EST output filename with downscaled values
  b1         =BONE       # opt: BUILD & EST slope filename
  b0         =BZERO      # opt: BUILD & EST intercept filename
  rsqr       =RSQR       # opt: BUILD regress. coeff. (r2) filename 
  resid      =RESID      # opt: BUILD residual (Yi - YHATi) filename
  ssxy       =SSXY       # opt: BUILD intermediate filename SSXY (sum of squares xy)
  ssxx       =SSXX       # opt: BUILD intermediate filename SSXX (sum of squares xx)
  ssyy       =SSYY       # opt: BUILD intermediate filename SSYY (sum of squares yy)
  d0         =1.0        # opt: BUILD wt value for 5 x 5 center cell      
  d1         =1.0        # opt: BUILD wt value at dist. = 1 cell from center
  d2         =1.0        # opt: BUILD wt value at dist. = 2 cells from center
  e0         =8.0        # opt: EST wt value for 5 x 5 center cell      
  e1         =2.0        # opt: EST wt value at dist. = 1 cell from center
  e2         =1.0        # opt: EST wt value at dist. = 2 cells from center
  rth        =0.08       # opt: EST minimum acceptable local rsqr from build
  #--------------------------------------------------------------------------------
  ######################################################################
END

# now, stuff all the above key=value pairs into a hash
%dscaleparams = ();
@allparams = split '\n', $allparams;
while(@allparams) {
    $_ = shift @allparams;
    s/#.*$//;  # delete comments
    if (($key, $value) = /(\S+)\s*=\s*(\S+)\s*$/) {
        if ($value eq "NULL") { $value = "" }
        $dscaleparams{$key} = $value;
#DEBUG print "key: $key value: $value \n";
    }
}

if ($#ARGV < 0 or $#ARGV == 0)  {   # try to process options
    if (($#ARGV < 0) or ($ARGV[0] =~ /(-h)|(--help)/))  {
        print<<END;
$0: Perl script to use r.mapcalc to downscale from low to high 
    resolution using 5 x 5 local regressions between variables
    like rainfall or temperature (dep. variable) and elevation 
    (indep. variable). 
    The program expects to be run twice: 
    Once at the lower resolution to build models using calc=BUILD,
    then at a higher resolution using calc=EST to calculate new 
    high resolution values. Example usage:

    GRASS >g.region res=0:02:30.00                             # set lower res
    GRASS >r.downscale.pl input=H_50 elev_lo=lodem calc=BUILD  # build model
    GRASS >g.region res=0:00:09.00                             # set higher res
    GRASS >r.downscale.pl elev_hi=DEM30m input=H_50 output=HI_50 calc=EST # est

    For help:    
    --help or -h        print this message and exit.
    --keys or -k        print the list of keys and exit.
    --defaults or -d    print the list of keys and defaults and exit.
   Other arguments are interpreted as parameter=value pairs

   Notes: For GRASS 5.0, no stats are scaled up: for grass 4.x find another
          version of this program that uses scaling. 
END
#print $#ARGV, "\n";
exit;        
    } elsif ($ARGV[0] =~ /(-k)|(--keys)/) {
        $count=0;
        for $key (sort keys %dscaleparams) {
            print "$key ";
            if (++$count > 5) {
                print "\n";
                $count=0;
            }
        }
        print "\n";
        exit;
        
    } elsif ($ARGV[0] =~ /(-)|(--defaults)/) {
        foreach $_ (split '\n', $allparams) {
            next if /##/;
            print "$_\n";
        }
        exit;
        
    }
}         # end of option processing


# now loop over command-line arguments
while($_ = shift @ARGV) {
    # parse them to key=value,
    ($key, $value) = /(.*)=(.*)/;

    # is the key defined?
    if (defined($dscaleparams{$key})) {
        if ($value eq "NULL") { $value = "" }     # "" can get dropped, so...
        $dscaleparams{$key}=$value;                # replace default with value.
    } else {
        # If not, it's probably a typo.
        die "$0: bad parameter '$key'\n";
    }
}

# DEBUG: print the hash
#for $key (sort keys %dscaleparams) {
#    print "$key\t=$dscaleparams{$key}\n";
#}
#exit;

# input parameters
my $x= $dscaleparams{'elev_lo'};
my $y= $dscaleparams{'input'};
my $hidem= $dscaleparams{'elev_hi'};
my $hiddtmp = "hiddtmp";
my $hidd= $dscaleparams{'output'};
my $calc= $dscaleparams{'calc'};
my $ssxy = $dscaleparams{'ssxy'};
my $ssxx = $dscaleparams{'ssxx'};
my $ssyy = $dscaleparams{'ssyy'};
my $rsqr = $dscaleparams{'rsqr'};
my $resid = $dscaleparams{'resid'};
my $b1 = $dscaleparams{'b1'}; 
my $b0 = $dscaleparams{'b0'};
# weights to use for wtd regression - 0 for middle cell 1 for dist.=1 cell, etc
# so for $w0=8, $w1=2 $w2=1, the mid cell gets 8X the weight of a cell 
#     2 cells distant, which would work for 1/D2, and $w0=3, $w1=1.5, 
#     and $w2=1 for 1/D
my $w0 = $dscaleparams{'d0'};
my $w1 = $dscaleparams{'d1'};
my $w2 = $dscaleparams{'d2'};
# and set of wts. for estimation using the model, to provide smoothing 
#  of results, these wts should possibly be different than model building,
# so for $e0=8, $e1=2 $e2=1, the mid cell gets 8X the weight of a 
#   cell 2 cells distant
my $e0 = $dscaleparams{'e0'};
my $e1 = $dscaleparams{'e1'};
my $e2 = $dscaleparams{'e2'};

my $rth = $dscaleparams{'rth'};

# constants - save a little typing
my $m12 = "[-2,-2]";
my $m11 = "[-2,-1]";
my $m10 = "[-2,0]";
my $m9 = "[-2,1]";
my $m8 = "[-2,2]";
my $m7 = "[-1,-2]";
my $m6 = "[-1,-1]";
my $m5 = "[-1,0]";
my $m4 = "[-1,1]";
my $m3 = "[-1,2]";
my $m2 = "[0,-2]";
my $m1 = "[0,-1]";
my $m0 = "[0,0]";
my $p1 = "[0,1]";
my $p2 = "[0,2]";
my $p3 = "[1,-2]";
my $p4 = "[1,-1]";
my $p5 = "[1,0]";
my $p6 = "[1,1]";
my $p7 = "[1,2]";
my $p8 = "[2,-2]";
my $p9 = "[2,-1]";
my $p10 = "[2,0]";
my $p11 = "[2,1]";
my $p12 = "[2,2]";

# use hash to specify different stored r.mapcalc procedures
%methods=(); 


# model building using weighted regression method on 5 x 5 local neighborhood
$methods{BUILD}=<<'EOL';
$ssxy = eval(xmean = ( if($x $m12<=0,$x $m0,$x $m12) * $w2 +  \ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 +  \ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 +  \ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 +  \ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 +  \ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 +  \ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 +  \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +  \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +  \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +  \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 +  \ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 +  \ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 +  \ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 +  \ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 +  \ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 +  \ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2\ eoln
       )/($w2 * 16 + $w1 * 8 + $w0),\ eoln
sm12 = if($x $m12 <=0,$x $m0 ,$x $m12) - xmean,\ eoln
sm11 = if($x $m11 <=0,$x $m0 ,$x $m11) - xmean,\ eoln
sm10 = if($x $m10 <=0,$x $m0 ,$x $m10) - xmean,\ eoln
sm9 = if($x $m9 <=0,$x $m0 ,$x $m9) - xmean,\ eoln
sm8 = if($x $m8 <=0,$x $m0 ,$x $m8) - xmean,\ eoln
sm7 = if($x $m7 <=0,$x $m0 ,$x $m7) - xmean,\ eoln
sm6 = if($x $m6 <=0,$x $m0 ,$x $m6) - xmean,\ eoln
sm5 = if($x $m5 <=0,$x $m0 ,$x $m5) - xmean,\ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4) - xmean,\ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3) - xmean,\ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean,\ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean,\ eoln
s0 = $x [0,0] - xmean,\ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1) - xmean,\ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2) - xmean,\ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3) - xmean,\ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4) - xmean,\ eoln
s5 = if($x $p5 <=0,$x $m0 ,$x $p5) - xmean,\ eoln
s6 = if($x $p6 <=0,$x $m0 ,$x $p6) - xmean,\ eoln
s7 = if($x $p7 <=0,$x $m0 ,$x $p7) - xmean,\ eoln
s8 = if($x $p8 <=0,$x $m0 ,$x $p8) - xmean,\ eoln
s9 = if($x $p9 <=0,$x $m0 ,$x $p9) - xmean,\ eoln
s10 = if($x $p10 <=0,$x $m0 ,$x $p10) - xmean,\ eoln
s11 = if($x $p11 <=0,$x $m0 ,$x $p11) - xmean,\ eoln
s12 = if($x $p12 <=0,$x $m0 ,$x $p12) - xmean,\ eoln
CR
       (sm12 * if($y $m12 <=0,$y $m0 ,$y $m12) * $w2 +\ eoln
       sm11 * if($y $m11 <=0,$y $m0 ,$y $m11) * $w2 +\ eoln
       sm10 * if($y $m10 <=0,$y $m0 ,$y $m10) * $w2 +\ eoln
       sm9 * if($y $m9 <=0,$y $m0 ,$y $m9) * $w2 +\ eoln
       sm8 * if($y $m8 <=0,$y $m0 ,$y $m8) * $w2 +\ eoln
       sm7 * if($y $m7 <=0,$y $m0 ,$y $m7) * $w2 +\ eoln
       sm6 * if($y $m6 <=0,$y $m0 ,$y $m6) * $w1 +\ eoln
       sm5 * if($y $m5 <=0,$y $m0 ,$y $m5) * $w1 +\ eoln
       sm4 * if($y $m4 <=0,$y $m0 ,$y $m4) * $w1 +\ eoln
       sm3 * if($y $m3 <=0,$y $m0 ,$y $m3) * $w2 +\ eoln
       sm2 * if($y $m2 <=0,$y $m0 ,$y $m2) * $w2 +\ eoln
       sm1 * if($y $m1 <=0,$y $m0 ,$y $m1) * $w1 +\ eoln
       s0 * ($y $m0) * $w0 +\ eoln
       s1 * if($y $p1 <=0,$y $m0 ,$y $p1) * $w1 +\ eoln
       s2 * if($y $p2 <=0,$y $m0 ,$y $p2) * $w2 +\ eoln
       s3 * if($y $p3 <=0,$y $m0 ,$y $p3) * $w2 +\ eoln
       s4 * if($y $p4 <=0,$y $m0 ,$y $p4) * $w1 +\ eoln
       s5 * if($y $p5 <=0,$y $m0 ,$y $p5) * $w1 +\ eoln
       s6 * if($y $p6 <=0,$y $m0 ,$y $p6) * $w1 +\ eoln
       s7 * if($y $p7 <=0,$y $m0 ,$y $p7) * $w2 +\ eoln
       s8 * if($y $p8 <=0,$y $m0 ,$y $p8) * $w2 +\ eoln
       s9 * if($y $p9 <=0,$y $m0 ,$y $p9) * $w2 +\ eoln
       s10 * if($y $p10 <=0,$y $m0 ,$y $p10) * $w2 +\ eoln
       s11 * if($y $p11 <=0,$y $m0 ,$y $p11) * $w2 +\ eoln
       s12 * if($y $p12 <=0,$y $m0 ,$y $p12) * $w2)) next
$ssxx = eval(\ eoln
xmean = (\ eoln
       if($x $m12<=0,$x $m0,$x $m12) * $w2 + \ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 + \ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 + \ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 + \ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 + \ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 + \ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 + \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 + \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 + \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 + \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 + \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 + \ eoln
       $x $m0 * $w0 + \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 + \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 + \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 + \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 + \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 + \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 + \ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 + \ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 + \ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 + \ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 + \ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 + \ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2\ eoln
       )/($w2 * 16 + $w1 * 8 + $w0), \ eoln
sm12 = if($x $m12 <=0,$x $m0 ,$x $m12 ) - xmean,\ eoln
sm11 = if($x $m11 <=0,$x $m0 ,$x $m11 ) - xmean,\ eoln
sm10 = if($x $m10 <=0,$x $m0 ,$x $m10 ) - xmean,\ eoln
sm9 = if($x $m9 <=0,$x $m0 ,$x $m9 ) - xmean,\ eoln
sm8 = if($x $m8 <=0,$x $m0 ,$x $m8 ) - xmean,\ eoln
sm7 = if($x $m7 <=0,$x $m0 ,$x $m7 ) - xmean,\ eoln
sm6 = if($x $m6 <=0,$x $m0 ,$x $m6 ) - xmean,\ eoln
sm5 = if($x $m5 <=0,$x $m0 ,$x $m5 ) - xmean,\ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4 ) - xmean,\ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3 ) - xmean,\ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean,\ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean,\ eoln
s0 = $x [0,0] - xmean,\ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1 ) - xmean,\ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2 ) - xmean,\ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3 ) - xmean,\ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4 ) - xmean,\ eoln
s5 = if($x $p5 <=0,$x $m0 ,$x $p5 ) - xmean,\ eoln
s6 = if($x $p6 <=0,$x $m0 ,$x $p6 ) - xmean,\ eoln
s7 = if($x $p7 <=0,$x $m0 ,$x $p7 ) - xmean,\ eoln
s8 = if($x $p8 <=0,$x $m0 ,$x $p8 ) - xmean,\ eoln
s9 = if($x $p9 <=0,$x $m0 ,$x $p9 ) - xmean,\ eoln
s10 = if($x $p10 <=0,$x $m0 ,$x $p10 ) - xmean,\ eoln
s11 = if($x $p11 <=0,$x $m0 ,$x $p11 ) - xmean,\ eoln
s12 = if($x $p12 <=0,$x $m0 ,$x $p12 ) - xmean,\ eoln
      (sm12 * sm12 * $w2 + sm11 * sm11 * $w2 + sm10 * sm10 * $w2 +\ eoln
      sm9 * sm9 * $w2 + sm8 * sm8 * $w2 + sm7 * sm7 * $w2 +\ eoln
      sm6 * sm6 * $w1 + sm5 * sm5 * $w1 + sm4 * sm4 * $w1 +\ eoln
      sm3 * sm3 * $w2 + sm2 * sm2 * $w2 + sm1 * sm1 * $w1 +\ eoln
      s0 * s0 * $w0 +\
      s1 * s1 * $w1 + s2 * s2 * $w2 + s3 * s3 * $w2 + s4 * s4 * $w1 +\ eoln
      s5 * s5 * $w1 + s6 * s6 * $w1 + s7 * s7 * $w2 + s8 * s8 * $w2 +\ eoln
      s9 * s9 * $w2 + s10 * s10 * $w2 + s11 * s11 * $w2 + s12 * s12 * $w2)\ eoln
) next
$b1 = (($ssxy) / ($ssxx)) next  
$b0 = eval( \ eoln
xmean = ( \ eoln
       if($x $m12<=0,$x $m0,$x $m12) * $w2 +\ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 +\ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 +\ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 +\ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 +\ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 +\ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 +\ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +\ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +\ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +\ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +\ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +\ eoln
       $x $m0 * $w0 +\ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +\ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +\ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +\ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +\ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +\ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 +\ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 +\ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 +\ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 +\ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 +\ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 +\ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2\ eoln
       )/($w2 * 16 + $w1 * 8 + $w0),\ eoln
ymean = (\ eoln
       if($y $m12 <=0,$y $m0 ,$y $m12 ) * $w2 +\ eoln
       if($y $m11 <=0,$y $m0 ,$y $m11 ) * $w2 +\ eoln
       if($y $m10 <=0,$y $m0 ,$y $m10 ) * $w2 +\ eoln
       if($y $m9 <=0,$y $m0 ,$y $m9 ) * $w2 +\ eoln
       if($y $m8 <=0,$y $m0 ,$y $m8 ) * $w2 +\ eoln
       if($y $m7 <=0,$y $m0 ,$y $m7 ) * $w2 +\ eoln
       if($y $m6 <=0,$y $m0 ,$y $m6 ) * $w1 +\ eoln
       if($y $m5 <=0,$y $m0 ,$y $m5 ) * $w1 +\ eoln
       if($y $m4 <=0,$y $m0 ,$y $m4 ) * $w1 +\ eoln
       if($y $m3 <=0,$y $m0 ,$y $m3 ) * $w2 +\ eoln
       if($y $m2 <=0,$y $m0 ,$y $m2 ) * $w2 +\ eoln
       if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 +\ eoln
       ($y $m0 ) * $w0 +\ eoln
       if($y $p1 <=0,$y $m0 ,$y $p1 ) * $w1 +\ eoln
       if($y $p2 <=0,$y $m0 ,$y $p2 ) * $w2 +\ eoln
       if($y $p3 <=0,$y $m0 ,$y $p3 ) * $w2 +\ eoln
       if($y $p4 <=0,$y $m0 ,$y $p4 ) * $w1 +\ eoln
       if($y $p5 <=0,$y $m0 ,$y $p5 ) * $w1 +\ eoln
       if($y $p6 <=0,$y $m0 ,$y $p6 ) * $w1 +\ eoln
       if($y $p7 <=0,$y $m0 ,$y $p7 ) * $w2 +\ eoln
       if($y $p8 <=0,$y $m0 ,$y $p8 ) * $w2 +\ eoln
       if($y $p9 <=0,$y $m0 ,$y $p9 ) * $w2 +\ eoln
       if($y $p10 <=0,$y $m0 ,$y $p10 ) * $w2 +\ eoln
       if($y $p11 <=0,$y $m0 ,$y $p11 ) * $w2 +\ eoln
       if($y $p12 <=0,$y $m0 ,$y $p12 ) * $w2 \ eoln
       )/($w2 * 16 + $w1 * 8 + $w0),\ eoln
ymean - ($b1 * xmean)  \ eoln
) next 
$ssyy = eval(\ eoln
ymean = (\ eoln
       if($y $m12 <=0,$y $m0 ,$y $m12 ) * $w2 +\ eoln
       if($y $m11 <=0,$y $m0 ,$y $m11 ) * $w2 +\ eoln
       if($y $m10 <=0,$y $m0 ,$y $m10 ) * $w2 +\ eoln
       if($y $m9 <=0,$y $m0 ,$y $m9 ) * $w2 +\ eoln
       if($y $m8 <=0,$y $m0 ,$y $m8 ) * $w2 +\ eoln
       if($y $m7 <=0,$y $m0 ,$y $m7 ) * $w2 +\ eoln
       if($y $m6 <=0,$y $m0 ,$y $m6 ) * $w1 +\ eoln
       if($y $m5 <=0,$y $m0 ,$y $m5 ) * $w1 +\ eoln
       if($y $m4 <=0,$y $m0 ,$y $m4 ) * $w1 +\ eoln
       if($y $m3 <=0,$y $m0 ,$y $m3 ) * $w2 +\ eoln
       if($y $m2 <=0,$y $m0 ,$y $m2 ) * $w2 +\ eoln
       if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 +\ eoln
       ($y $m0 ) * $w0 +\ eoln
       if($y $p1 <=0,$y $m0 ,$y $p1 ) * $w1 +\ eoln
       if($y $p2 <=0,$y $m0 ,$y $p2 ) * $w2 +\ eoln
       if($y $p3 <=0,$y $m0 ,$y $p3 ) * $w2 +\ eoln
       if($y $p4 <=0,$y $m0 ,$y $p4 ) * $w1 +\ eoln
       if($y $p5 <=0,$y $m0 ,$y $p5 ) * $w1 +\ eoln
       if($y $p6 <=0,$y $m0 ,$y $p6 ) * $w1 +\ eoln
       if($y $p7 <=0,$y $m0 ,$y $p7 ) * $w2 +\ eoln
       if($y $p8 <=0,$y $m0 ,$y $p8 ) * $w2 +\ eoln
       if($y $p9 <=0,$y $m0 ,$y $p9 ) * $w2 +\ eoln
       if($y $p10 <=0,$y $m0 ,$y $p10 ) * $w2 +\ eoln
       if($y $p11 <=0,$y $m0 ,$y $p11 ) * $w2 +\ eoln
       if($y $p12 <=0,$y $m0 ,$y $p12 ) * $w2\ eoln
       )/($w2 * 16 + $w1 * 8 + $w0),\ eoln
sm12 = if($y $m12 <=0,$y $m0 ,$y $m12 ) - ymean,\ eoln
sm11 = if($y $m11 <=0,$y $m0 ,$y $m11 ) - ymean,\ eoln
sm10 = if($y $m10 <=0,$y $m0 ,$y $m10 ) - ymean,\ eoln
sm9 = if($y $m9 <=0,$y $m0 ,$y $m9 ) - ymean,\ eoln
sm8 = if($y $m8 <=0,$y $m0 ,$y $m8 ) - ymean,\ eoln
sm7 = if($y $m7 <=0,$y $m0 ,$y $m7 ) - ymean,\ eoln
sm6 = if($y $m6 <=0,$y $m0 ,$y $m6 ) - ymean,\ eoln
sm5 = if($y $m5 <=0,$y $m0 ,$y $m5 ) - ymean,\ eoln
sm4 = if($y $m4 <=0,$y $m0 ,$y $m4 ) - ymean,\ eoln
sm3 = if($y $m3 <=0,$y $m0 ,$y $m3 ) - ymean,\ eoln
sm2 = if($y $m2 <=0,$y $m0 ,$y $m2 ) - ymean,\ eoln
sm1 = if($y $m1 <=0,$y $m0 ,$y $m1 ) - ymean,\ eoln
s0 = if($y $m0 <=0,$y $m0 ,$y $m0 ) - ymean,\ eoln
s1 = if($y $p1 <=0,$y $m0 ,$y $p1 ) - ymean,\ eoln
s2 = if($y $p2 <=0,$y $m0 ,$y $p2 ) - ymean,\ eoln
s3 = if($y $p3 <=0,$y $m0 ,$y $p3 ) - ymean,\ eoln
s4 = if($y $p4 <=0,$y $m0 ,$y $p4 ) - ymean,\ eoln
s5 = if($y $p5 <=0,$y $m0 ,$y $p5 ) - ymean,\ eoln
s6 = if($y $p6 <=0,$y $m0 ,$y $p6 ) - ymean,\ eoln
s7 = if($y $p7 <=0,$y $m0 ,$y $p7 ) - ymean,\ eoln
s8 = if($y $p8 <=0,$y $m0 ,$y $p8 ) - ymean,\ eoln
s9 = if($y $p9 <=0,$y $m0 ,$y $p9 ) - ymean,\ eoln
s10 = if($y $p10 <=0,$y $m0 ,$y $p10 ) - ymean,\ eoln
s11 = if($y $p11 <=0,$y $m0 ,$y $p11 ) - ymean,\ eoln
s12 = if($y $p12 <=0,$y $m0 ,$y $p12 ) - ymean,\ eoln
(sm12 * sm12 * $w2 + sm11 * sm11 * $w2 + sm10 * sm10 * $w2 +\ eoln
sm9 * sm9 * $w2 + sm8 * sm8 * $w2 + sm7 * sm7 * $w2 +\ eoln
sm6 * sm6 * $w1 + sm5 * sm5 * $w1 + sm4 * sm4 * $w1 +\ eoln
sm3 * sm3 * $w2 + sm2 * sm2 * $w2 + sm1 * sm1 * $w1 +\ eoln
s0 * s0 * $w0 +\ eoln
s1 * s1 * $w1 + s2 * s2 * $w2 + s3 * s3 * $w2 + s4 * s4 * $w1 +\ eoln
s5 * s5 * $w1 + s6 * s6 * $w1 + s7 * s7 * $w2 + s8 * s8 * $w2 +\ eoln
s9 * s9 * $w2 + s10 * s10 * $w2 + s11 * s11 * $w2 + s12 * s12 * $w2)) next
$rsqr = eval (rval = $ssxy/ sqrt($ssxx * $ssyy), \ eoln
(rval * rval)) next
$resid = $y - ($x * $b1 + $b0) next
ANOM = abs($resid)/$rsqr \ eoln
EOL

#calculate estimated DDs from regression params - simple
$methods{EST3}=<<'EOL';
$hidd = $hidem * $b1 + $b0 \ eoln
EOL

#calculate estimated DDs from regression params - weighted avg neighbor 5x5 only
# weighting approximate discrete version of 1/distance-squared 8-2-1 for distance
# 1-2-3
$methods{EST9}=<<'EOL';
$hidd = (if($rsqr $m6 <=$rth,$y $m6, \ eoln
              $b0 $m6 + ($hidem * $b1 $m6)) * $e1 + \ eoln
        if($rsqr $m5 <=$rth,$y $m5, \ eoln
              $b0 $m5 + ($hidem * $b1 $m5)) * $e1 + \ eoln
        if($rsqr $m4 <=$rth,$y $m4, \ eoln
              $b0 $m4 + ($hidem * $b1 $m4)) * $e1 + \ eoln
        if($rsqr $m1 <=$rth,$y $m1, \ eoln
              $b0 $m1 + ($hidem * $b1 $m1)) * $e1 + \ eoln
        if($rsqr $m0 <=$rth,$y $m0, \ eoln
              $b0 $m0 + ($hidem * $b1 $m0)) * $e0 + \ eoln
        if($rsqr $p1 <=$rth,$y $p1, \ eoln
              $b0 $p1 + ($hidem * $b1 $p1)) * $e1 + \ eoln
        if($rsqr $p4 <=$rth,$y $p4, \ eoln
              $b0 $p4 + ($hidem * $b1 $p4)) * $e1 + \ eoln
        if($rsqr $p5 <=$rth,$y $p5, \ eoln
              $b0 $p5 + ($hidem * $b1 $p5)) * $e1 + \ eoln
        if($rsqr $p6 <=$rth,$y $p6, \ eoln
              $b0 $p6 + ($hidem * $b1 $p6)) * $e1) \ eoln
        /($e1 * 8 + $e0) \ eoln
EOL
$methods{EST}=<<'EOL';
$hidd = (if($rsqr $m12 <=$rth,$y $m12, \ eoln
              $b0 $m12 + ($hidem * $b1 $m12)) * $e2 + \ eoln
        if($rsqr $m11 <=$rth,$y $m11, \ eoln
              $b0 $m11 + ($hidem * $b1 $m11)) * $e2 + \ eoln
        if($rsqr $m10 <=$rth,$y $m10, \ eoln
              $b0 $m10 + ($hidem * $b1 $m10)) * $e2 + \ eoln
        if($rsqr $m9 <=$rth,$y $m9, \ eoln
              $b0 $m9 + ($hidem * $b1 $m9)) * $e2 + \ eoln
        if($rsqr $m8 <=$rth,$y $m8, \ eoln
              $b0 $m8 + ($hidem * $b1 $m8)) * $e2 + \ eoln
        if($rsqr $m7 <=$rth,$y $m7, \ eoln
              $b0 $m7 + ($hidem * $b1 $m7)) * $e2 + \ eoln
        if($rsqr $m6 <=$rth,$y $m6, \ eoln
              $b0 $m6 + ($hidem * $b1 $m6)) * $e1 + \ eoln
        if($rsqr $m5 <=$rth,$y $m5, \ eoln
              $b0 $m5 + ($hidem * $b1 $m5)) * $e1 + \ eoln
        if($rsqr $m4 <=$rth,$y $m4, \ eoln
              $b0 $m4 + ($hidem * $b1 $m4)) * $e1 + \ eoln
        if($rsqr $m3 <=$rth,$y $m3, \ eoln
              $b0 $m3 + ($hidem * $b1 $m3)) * $e2 + \ eoln
        if($rsqr $m2 <=$rth,$y $m2, \ eoln
              $b0 $m2 + ($hidem * $b1 $m2)) * $e2 + \ eoln
        if($rsqr $m1 <=$rth,$y $m1, \ eoln
              $b0 $m1 + ($hidem * $b1 $m1)) * $e1 + \ eoln
        if($rsqr $m0 <=$rth,$y $m0, \ eoln
              $b0 $m0 + ($hidem * $b1 $m0)) * $e0 + \ eoln
        if($rsqr $p1 <=$rth,$y $p1, \ eoln
              $b0 $p1 + ($hidem * $b1 $p1)) * $e1 + \ eoln
        if($rsqr $p2 <=$rth,$y $p2, \ eoln
              $b0 $p2 + ($hidem * $b1 $p2)) * $e2 + \ eoln
        if($rsqr $p3 <=$rth,$y $p3, \ eoln
              $b0 $p3 + ($hidem * $b1 $p3)) * $e2 + \ eoln
        if($rsqr $p4 <=$rth,$y $p4, \ eoln
              $b0 $p4 + ($hidem * $b1 $p4)) * $e1 + \ eoln
        if($rsqr $p5 <=$rth,$y $p5, \ eoln
              $b0 $p5 + ($hidem * $b1 $p5)) * $e1 + \ eoln
        if($rsqr $p6 <=$rth,$y $p6, \ eoln
              $b0 $p6 + ($hidem * $b1 $p6)) * $e1 + \ eoln
        if($rsqr $p7 <=$rth,$y $p7, \ eoln
              $b0 $p7 + ($hidem * $b1 $p7)) * $e2 + \ eoln
        if($rsqr $p8 <=$rth,$y $p8, \ eoln
              $b0 $p8 + ($hidem * $b1 $p8)) * $e2 + \ eoln
        if($rsqr $p9 <=$rth,$y $p9, \ eoln
              $b0 $p9 + ($hidem * $b1 $p9)) * $e2 + \ eoln
        if($rsqr $p10 <=$rth,$y $p10, \ eoln
              $b0 $p10 + ($hidem * $b1 $p10)) * $e2 + \ eoln
        if($rsqr $p11 <=$rth,$y $p11, \ eoln
              $b0 $p11 + ($hidem * $b1 $p11)) * $e2 + \ eoln
        if($rsqr $p12 <=$rth,$y $p12, \ eoln
              $b0 $p12 + ($hidem * $b1 $p12)) * $e2) \ eoln
        /($e2 * 16 + $e1 * 8 + $e0) \ eoln
EOL
$methods{ESTold}=<<'EOL';
$hidd = (if($b0 $m12 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m12 + ($hidem * $b1 $m12)) * $e2 + \ eoln
        if($b0 $m11 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m11 + ($hidem * $b1 $m11)) * $e2 + \ eoln
        if($b0 $m10 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m10 + ($hidem * $b1 $m10)) * $e2 + \ eoln
        if($b0 $m9 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m9 + ($hidem * $b1 $m9)) * $e2 + \ eoln
        if($b0 $m8 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m8 + ($hidem * $b1 $m8)) * $e2 + \ eoln
        if($b0 $m7 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m7 + ($hidem * $b1 $m7)) * $e2 + \ eoln
        if($b0 $m6 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m6 + ($hidem * $b1 $m6)) * $e1 + \ eoln
        if($b0 $m5 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m5 + ($hidem * $b1 $m5)) * $e1 + \ eoln
        if($b0 $m4 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m4 + ($hidem * $b1 $m4)) * $e1 + \ eoln
        if($b0 $m3 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m3 + ($hidem * $b1 $m3)) * $e2 + \ eoln
        if($b0 $m2 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m2 + ($hidem * $b1 $m2)) * $e2 + \ eoln
        if($b0 $m1 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $m1 + ($hidem * $b1 $m1)) * $e1 + \ eoln
              ($b0 $m0 + ($hidem * $b1 $m0)) * $e0 + \ eoln
        if($b0 $p1 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p1 + ($hidem * $b1 $p1)) * $e1 + \ eoln
        if($b0 $p2 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p2 + ($hidem * $b1 $p2)) * $e2 + \ eoln
        if($b0 $p3 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p3 + ($hidem * $b1 $p3)) * $e2 + \ eoln
        if($b0 $p4 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p4 + ($hidem * $b1 $p4)) * $e1 + \ eoln
        if($b0 $p5 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p5 + ($hidem * $b1 $p5)) * $e1 + \ eoln
        if($b0 $p6 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p6 + ($hidem * $b1 $p6)) * $e1 + \ eoln
        if($b0 $p7 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p7 + ($hidem * $b1 $p7)) * $e2 + \ eoln
        if($b0 $p8 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p8 + ($hidem * $b1 $p8)) * $e2 + \ eoln
        if($b0 $p9 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p9 + ($hidem * $b1 $p9)) * $e2 + \ eoln
        if($b0 $p10 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p10 + ($hidem * $b1 $p10)) * $e2 + \ eoln
        if($b0 $p11<=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p11 + ($hidem * $b1 $p11)) * $e2 + \ eoln
        if($b0 $p12 <=0,$b0 $m0 + ($hidem * $b1 $m0), \ eoln
              $b0 $p12 + ($hidem * $b1 $p12)) * $e2) \ eoln
        /($e2 * 16 + $e1 * 8 + $e0) \ eoln
EOL

#  test simple avg
$methods{R}='AVG = ($x [-1,-1] + $x [0,0] + $x [1,1])/3';
# test for eval and eoln rigamorall
$methods{R1}=<<'EOL';
AVG = eval( \ eoln
t1=$x [-1,-1], \ eoln
t2=$x [0,0], \ eoln
t3=$x [1,1], \ eoln
(t1+t2+t3)/3)
EOL
# ratio method DDhires = ELEVhires x DDlores / ELEVlores
$methods{R2}=<<'EOL';
#RATIO = $y / $x  eoln 
#OTHER = $x / $y
$ssxy = eval( \ eoln
xmean = (if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
            $x $m0 * $w0 + \ eoln
         if($x $p1<=0,$x $m0,$x $p1) * $w1 )/3, \ eoln
CR
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean, \ eoln
s0 = $x [0,0] - xmean, \ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1 ) - xmean, \ eoln
CR
(sm1 * if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 + \ eoln
s0 * ($y $m0 ) * $w0 + \ eoln
EOL

# ratio method DDhires = ELEVhires x DDlores / ELEVlores
$methods{R3}=<<'EOL';
$ssxy = eval(xmean = ( 
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 \ eoln
       )/($w2 * 4 + $w1 * 4 + $w0),\ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4) - xmean,\ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3) - xmean,\ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean,\ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean,\ eoln
s0 = $x [0,0] - xmean,\ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1) - xmean,\ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2) - xmean,\ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3) - xmean,\ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4) - xmean,\ eoln
CR
       (sm4 * if($y $m4 <=0,$y $m0 ,$y $m4) * $w1 +\ eoln
       sm3 * if($y $m3 <=0,$y $m0 ,$y $m3) * $w2 +\ eoln
       sm2 * if($y $m2 <=0,$y $m0 ,$y $m2) * $w2 +\ eoln
       sm1 * if($y $m1 <=0,$y $m0 ,$y $m1) * $w1 +\ eoln
       s0 * ($y $m0) * $w0 +\ eoln
       s1 * if($y $p1 <=0,$y $m0 ,$y $p1) * $w1 +\ eoln
       s2 * if($y $p2 <=0,$y $m0 ,$y $p2) * $w2 +\ eoln
       s3 * if($y $p3 <=0,$y $m0 ,$y $p3) * $w2 + \ eoln
       s4 * if($y $p4 <=0,$y $m0 ,$y $p4) * $w1)) \ eoln
EOL
#RATIO = $y / $x  eoln 
#OTHER = $x / $y

# ratio method DDhires = ELEVhires x DDlores / ELEVlores
$methods{R4}=<<'EOL';
$ssxy = eval(xmean = (if($x $m6<=0,$x $m0,$x $m6) * $w1 +  \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +  \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +  \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +  \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1\ eoln
       )/($w2 * 4 + $w1 * 8 + $w0),\ eoln
sm6 = if($x $m6 <=0,$x $m0 ,$x $m6) - xmean,\ eoln
sm5 = if($x $m5 <=0,$x $m0 ,$x $m5) - xmean,\ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4) - xmean,\ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3) - xmean,\ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean,\ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean,\ eoln
s0 = $x [0,0] - xmean,\ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1) - xmean,\ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2) - xmean,\ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3) - xmean,\ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4) - xmean,\ eoln
s5 = if($x $p5 <=0,$x $m0 ,$x $p5) - xmean,\ eoln
s6 = if($x $p6 <=0,$x $m0 ,$x $p6) - xmean,\ eoln
CR
       (sm6 * if($y $m6 <=0,$y $m0 ,$y $m6) * $w1 +\ eoln
       sm5 * if($y $m5 <=0,$y $m0 ,$y $m5) * $w1 +\ eoln
       sm4 * if($y $m4 <=0,$y $m0 ,$y $m4) * $w1 +\ eoln
       sm3 * if($y $m3 <=0,$y $m0 ,$y $m3) * $w2 +\ eoln
       sm2 * if($y $m2 <=0,$y $m0 ,$y $m2) * $w2 +\ eoln
       sm1 * if($y $m1 <=0,$y $m0 ,$y $m1) * $w1 +\ eoln
       s0 * ($y $m0) * $w0 +\ eoln
       s1 * if($y $p1 <=0,$y $m0 ,$y $p1) * $w1 +\ eoln
       s2 * if($y $p2 <=0,$y $m0 ,$y $p2) * $w2 +\ eoln
       s3 * if($y $p3 <=0,$y $m0 ,$y $p3) * $w2 +\ eoln
       s4 * if($y $p4 <=0,$y $m0 ,$y $p4) * $w1 +\ eoln
       s5 * if($y $p5 <=0,$y $m0 ,$y $p5) * $w1 +\ eoln
       s6 * if($y $p6 <=0,$y $m0 ,$y $p6) * $w1)) \ eoln
EOL


# clean up r.mapcalc input template: $methods{$calc}
$methods{$calc} =~ s/#.*?$//mg;
$methods{$calc} =~ s/\n//mg;
$methods{$calc} =~ s/\s\s//mg;
$methods{$calc} =~ s/CR//mg;
$methods{$calc} =~ s/\s+eoln/ eoln\n/mg;
$methods{$calc} =~ s/next/ \n/mg;
$methods{$calc} =~ s/\\//mg;

# ----- set up r.mapcalc pipe -----
# sed required because perl cannot allow final slash marks expected by r.mapcalc
open RMAPCALC, "| sed -e 's|eoln|\\\\|g' | r.mapcalc "
  or die "Can't launch r.mapcalc";
#print "prior to eval\n $methods{$calc}\n";
 # trickery to interpolate the parameters into the command string
    # i.e. to replace "$var" in with the current value of $var in the string
   $calcstring =  eval '"' . "$methods{$calc}" . '\n"';

# use for debug:
#print "post eval\n";
#print "$calcstring\n";

# ----- actual call to r.mapcalc passing in the calc method -----
print RMAPCALC $calcstring;
close RMAPCALC;
# end of r.downscale.pl