set DATES ordered;

param hi {DATES};
param avg {DATES};
param lo {DATES};

param pi := 4*atan(1);

var a {j in 0..3};

var dev {DATES} >= 0, := 1;

minimize sumdev: sum {d in DATES} dev[d];

subject to def_pos_dev {d in DATES}:
    a[0] + a[1]*ord(d,DATES) + a[2]*cos( 2*pi*ord(d,DATES)/365.25)
                             + a[3]*sin( 2*pi*ord(d,DATES)/365.25)
	- avg[d] 
    <= dev[d];

subject to def_neg_dev {d in DATES}:
    -dev[d] <= 
    a[0] + a[1]*ord(d,DATES) + a[2]*cos( 2*pi*ord(d,DATES)/365.25)
                             + a[3]*sin( 2*pi*ord(d,DATES)/365.25)
	- avg[d];

data;
set DATES := include "dates.txt";
param: hi avg lo := include "WXDailyHistory.txt";

let a[0] := 60;
let a[1] := 0;
let a[2] := 20;
let a[3] := 20;

solve;

display a;

display a[1]*365.25;

#printf {d in DATES}: "%s %d %f %f \n", d, ord(d,DATES), avg[d], 
#    a[0] + a[1]*ord(d,DATES) + a[2]*cos( 2*pi*ord(d,DATES)/365.25)
#                             + a[3]*sin( 2*pi*ord(d,DATES)/365.25)
#	> "output";

#printf {d in DATES}: "L\t%f\t%f\t%f\n", 2005+(ord(d,DATES)-136)/365.25, avg[d], 
#    a[0] + a[1]*ord(d,DATES) + a[2]*cos( 2*pi*ord(d,DATES)/365.25)
#                             + a[3]*sin( 2*pi*ord(d,DATES)/365.25)
#	> "svg_output";

#printf {d in DATES}: "%f\t%f\t%f\n", 2005+(ord(d,DATES)-136)/365.25, avg[d], 
#    a[0] + a[1]*ord(d,DATES) + a[2]*cos( 2*pi*ord(d,DATES)/365.25)
#                             + a[3]*sin( 2*pi*ord(d,DATES)/365.25)
#	> "mat_output";