program analysis;
{
  This program calculates all available station derivatives each
  year and their average to estimate the global derivative. First
  the program reads in the lines of GYA.txt, the yearly averages for
  all station. It sorts these in order of station and year. It turns
  out that our data file was pre-sorted.  
}

const
  max_num_measurements=1000000;

type
  measurement_type=record
    id:longint;{station identifier}
    year:integer;{year of the measurement}
    temperature:real;{year-long average}
    derivative_valid:boolean;{valid if last year's measurement exists}
    derivative:real;{change from last year}
  end;
  
var
  data:text;{the data file}
  measurements:array [1..max_num_measurements] of measurement_type;
  measurement_num,num_measurements,y,counter:integer;
  i:longint;
  t,sum:real;
  sorted:boolean;
  
{
  We use this swap procedure to sort our list of stations. We
  use a "bubble sort" because it keeps the code simple.
}
procedure swap(p,q:integer);
var
  m:measurement_type;
begin
  m:=measurements[p];
  measurements[p]:=measurements[q];
  measurements[q]:=m;
end;

{
  This is the main program.
}
begin
{
  Start by reading in the Global Yearly Averages from our GYA.txt file,
  which we have formatted for our convenience, but which contains raw 
  data from NCDC. Each line in the file provides us with a station ID,
  a year, and a temperature.
}
  writeln('Read GYA.txt...');
  reset(data,'GYA.txt');
  num_measurements:=0;
  while not eof(data) do begin
    inc(num_measurements);
    with measurements[num_measurements] do 
      readln(data,id,year,temperature);
  end;
  close(data);
{
	We sort the measurements by increasing station ID and increasing 
	year.
}
  writeln('Sort measurements by station and year...');
  sorted:=false;
  while not sorted do begin
    sorted:=true;
    for measurement_num:=1 to num_measurements-1 do begin
      if (measurements[measurement_num].id > measurements[measurement_num+1].id)
                    or
        ((measurements[measurement_num].id = measurements[measurement_num+1].id)
                    and
        (measurements[measurement_num].year > measurements[measurement_num+1].year) )
                  then begin
          swap(measurement_num,measurement_num+1);
          sorted:=false;
      end;
    end;
  end;
{
  Each measurement now gets a derivative, provided that the previous measurement
  is from the same station in the previous year.
}
  writeln('Calculate and mark all valid derivatives...');  
  measurements[1].derivative_valid:=false;
  for measurement_num:=2 to num_measurements do begin
    if ((measurements[measurement_num].id = measurements[measurement_num-1].id)
                    and
    (measurements[measurement_num].year - 1 = measurements[measurement_num-1].year) )
                  then begin
      measurements[measurement_num].derivative_valid:=true;
      measurements[measurement_num].derivative:=
        measurements[measurement_num].temperature
        - measurements[measurement_num-1].temperature;
    end else begin
      measurements[measurement_num].derivative_valid:=false;
      measurements[measurement_num].derivative:=-999;
    end;
{    with measurements[measurement_num] do writeln(id:1,' ',year:1,' ',derivative:10:6);}
  end;

{
  For each year in our time range, calculate the average derivative of temperature
  and print it to the screen.
}
  writeln('Calculate average derivative for each year...');
  for y:=1800 to 2010 do begin
    sum:=0;
    counter:=0;
    for measurement_num:=1 to num_measurements do
      with measurements[measurement_num] do
        if (year=y) and derivative_valid then begin
          sum:=sum+derivative;
          inc(counter);
        end;
    writeln(y:1,' ',sum/counter:10:6,' ',counter:1);
  end;
{
  What you see now in the console is the derivatives with year. We pasted these into
  Excel and integrated them. You will find our Climate.xls spreadsheet that performs the
  integration in our Climate.zip archive.
  
  http://www.hashemifamily.com/Kevan/Climate/Climate.zip
  
  As you can see, the integration proceeds by adding the previous year's integral to
  the current year's average derivative. We perform no smoothing.
}
end.