## Octave plots for Mike's NYT Trick
## Jean S, 2014

clear all; close all;

graphics_toolkit('fltk');

## Download Mann's filter and MBH98 data 
if !exist("lowmann.m","file") 
  urlwrite("http://www.climateaudit.info/data/uc/lowmann.m",
           "lowmann.m");  
end
if !exist("nhmean.txt","file") 
  urlwrite("http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt",
           "nhmean.txt");  
end

## YEAR    NH RECON    NH RAW    +1 sigma   -1sigma   +2 sigma    -2 sigma   [file labels wrong ...]
MBH98=dlmread("nhmean.txt","",1,0);

I98_recon=MBH98(:,1)>=1400 & MBH98(:,1)<=1980;
I98_instru=MBH98(:,1)>=1902 & MBH98(:,1)<=1995;
I98_pre_instru=MBH98(:,1)>=1400 & MBH98(:,1)<=1901;
I98_i_splice=MBH98(:,1)>=1981 & MBH98(:,1)<=1995;
I98_smooth=MBH98(:,1)>=1425 & MBH98(:,1)<=1955;
I98_smooth_wt=MBH98(:,1)>=1425 & MBH98(:,1)<=1970;

smooth98=lowmann(MBH98(I98_recon,2),2/50);
smooth_w_trick98=lowmann([MBH98(I98_recon,2);MBH98(I98_i_splice,3)],2/50);
## Cut 25 samples from the both ends
smooth98=smooth98(26:end-25);
smooth_w_trick98=smooth_w_trick98(26:end-25);

## Emulation of the MBH98 Fig 5b w/o and w/ the trick
figure(1);
subplot(1,2,1); hold on;
plot(MBH98(I98_recon,1),MBH98(I98_recon,6),"linestyle",":","linewidth",1,'k');
plot(MBH98(I98_recon,1),MBH98(I98_recon,7),"linestyle",":","linewidth",1,'k');
plot(MBH98(I98_recon,1),MBH98(I98_recon,2),"linestyle","-","linewidth",1,'k');
plot(MBH98(I98_instru,1),MBH98(I98_instru,3),"linestyle","--","linewidth",1,'k');
plot(MBH98(I98_smooth,1),smooth98,"linestyle","--","linewidth",4,'r');
title("MBH98 Fig 5b without the trick");
box on;
set(gca,"xtick",1400:50:2000,"ytick",-0.8:0.2:0.6);
axis([1850 1997 -0.8 0.7]);
subplot(1,2,2); hold on;
plot(MBH98(I98_recon,1),MBH98(I98_recon,6),"linestyle",":","linewidth",1,'k');
plot(MBH98(I98_recon,1),MBH98(I98_recon,7),"linestyle",":","linewidth",1,'k');
plot(MBH98(I98_recon,1),MBH98(I98_recon,2),"linestyle","-","linewidth",1,'k');
plot(MBH98(I98_instru,1),MBH98(I98_instru,3),"linestyle","--","linewidth",1,'k');
plot(MBH98(I98_smooth_wt,1),smooth_w_trick98,"linestyle","--","linewidth",4,'r');
title("MBH98 Fig 5b with the trick");
box on;
set(gca,"xtick",1400:50:2000,"ytick",-0.8:0.2:0.6);
axis([1850 1997 -0.8 0.7])


## Download MBH99 data 
if !exist("nhem-raw.dat","file") 
  urlwrite("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann1999/recons/nhem-raw.dat",
           "nhem-raw.dat");  
end
if !exist("nhem-recon.dat","file") 
  urlwrite("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann1999/recons/nhem-recon.dat",
           "nhem-recon.dat");  
end
if !exist("nhem-standerr-labeled.dat","file") 
  urlwrite("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/mann1999/recons/nhem-standerr-labeled.dat",
           "nhem-standerr-labeled.dat");  
end

## YEAR    NH RECON
MBH99=dlmread("nhem-recon.dat","",0,0);
 ## YEAR   1 SIGMA    2 SIGMA   IGNORE THESE COLUMNS
MBH99_CIs=dlmread("nhem-standerr-labeled.dat","",1,0);
## YEAR    NH INSTRU
MBH99_raw=dlmread("nhem-raw.dat","",1,0);

I99_i_splice=MBH99_raw(:,1)>=1981 & MBH99_raw(:,1)<=1997;
I98_instru=MBH99_raw(:,1)>=1902 & MBH99_raw(:,1)<=1997;
I99_CIs=MBH99(:,1)>=1000 & MBH99(:,1)<=1901;
I99_smooth=MBH99(:,1)>=1000 & MBH99(:,1)<=1960;
I99_smooth_wt=MBH99(:,1)>=1000 & MBH99(:,1)<=1977;
I99_mannkovitch=MBH99(:,1)>=1000 & MBH99(:,1)<=1850;

smooth99=lowmann(MBH99(:,2),2/40);
smooth_w_trick99=lowmann([MBH99(:,2);MBH99_raw(I99_i_splice,2)],2/40);
## Cut 20 samples from the end
smooth99=smooth99(1:end-20);
smooth_w_trick99=smooth_w_trick99(1:end-20);

## Emulation of the MBH99 Fig 3a w/o and w/ the trick
figure(2); 
subplot(1,2,1); hold on;
fill([MBH99_CIs(:,1);flipud(MBH99_CIs(:,1))],...
        [MBH99(I99_CIs,2)-MBH99_CIs(:,3);flipud(MBH99(I99_CIs,2)+MBH99_CIs(:,3))],...
        "y","EdgeColor","None");
plot(MBH99(I99_mannkovitch,1),MBH99(I99_mannkovitch,2)-detrend(MBH99(I99_mannkovitch,2)),...
        "linestyle","--","linewidth",1,"k");
plot([1000,2000],[0,0],"linestyle","--","k");
plot([1000,2000],[MBH99_raw(end,2),MBH99_raw(end,2)],"linestyle",":","k");
plot(MBH99_raw(:,1),MBH99_raw(:,2),"linestyle","-","linewidth",1,"r");
plot(MBH99_raw(end,1),MBH99_raw(end,2),"markersize",18,"kx");
plot(MBH99(:,1),MBH99(:,2),"linestyle","-","linewidth",1,"k");
plot(MBH99(I99_smooth,1),smooth99,"linestyle","-","linewidth",3,"color",[0.5 0 0.7]);
title("MBH99 Fig 3a without the trick");
box on;
set(gca,"xtick",1000:100:2000,"ytick",-1.0:0.25:1.0);
axis([1850 2000 -1 1]);
subplot(1,2,2); hold on;
fill([MBH99_CIs(:,1);flipud(MBH99_CIs(:,1))],...
        [MBH99(I99_CIs,2)-MBH99_CIs(:,3);flipud(MBH99(I99_CIs,2)+MBH99_CIs(:,3))],...
        "y","EdgeColor","None");
plot(MBH99(I99_mannkovitch,1),MBH99(I99_mannkovitch,2)-detrend(MBH99(I99_mannkovitch,2)),...
        "linestyle","--","linewidth",1,'k');
plot([1000,2000],[0,0],"linestyle","--","k");
plot([1000,2000],[MBH99_raw(end,2),MBH99_raw(end,2)],"linestyle",":","k");
plot(MBH99_raw(:,1),MBH99_raw(:,2),"linestyle","-","linewidth",1,"r");
plot(MBH99_raw(end,1),MBH99_raw(end,2),"markersize",18,"kx");
plot(MBH99(:,1),MBH99(:,2),"linestyle","-","linewidth",1,"k");
plot(MBH99(I99_smooth_wt,1),smooth_w_trick99,"linestyle","-","linewidth",3,"color",[0.5 0 0.7]);
title("MBH99 Fig 3a with the trick");
box on;
set(gca,"xtick",1000:100:2000,"ytick",-1.0:0.25:1.0);
axis([1850 2000 -1 1]); 


## Emulation of the MBH98 Press Release Figure

figure(3); hold on;
fill([MBH98(I98_recon,1);flipud(MBH98(I98_recon,1))],...
        [MBH98(I98_recon,6);flipud(MBH98(I98_recon,7))],...
        "c","EdgeColor","None");
plot(MBH98(I98_recon,1),MBH98(I98_recon,2),"linestyle","-","linewidth",1,"k");
plot(MBH98(I98_smooth_wt,1),smooth_w_trick98,"linestyle","-","linewidth",4,"k");
plot(MBH99_raw(I99_i_splice,1),MBH99_raw(I99_i_splice,2),"linestyle","-","linewidth",1,"r");
title("MBH98 Press Release Figure 2")
axis([1400 2000 -0.8 0.7])
set(gca,"xtick",1400:200:2000,"ytick",-0.6:0.3:0.6);


## Emulation of the NYT Figure w/o and w/ the splicing 

figure(4); 
subplot(2,1,1); hold on;
plot(MBH98(I98_recon,1),MBH98(I98_recon,2)*1.8,"linestyle","-","linewidth",1,'k');
plot([1400,2000],[0,0],"linestyle","-","k");
plot([1400,2000],[0.2,0.2],"linestyle",":","k");
plot([1400,2000],[0.4,0.4],"linestyle",":","k");
plot([1400,2000],[-0.2,-0.2],"linestyle",":","k");
plot([1400,2000],[-0.4,-0.4],"linestyle",":","k");
plot([1400,2000],[-0.6,-0.6],"linestyle",":","k");
plot([1400,2000],[-0.8,-0.8],"linestyle",":","k");
plot([1900,2000],[0.6,0.6],"linestyle",":","k");
plot([1925,2000],[0.8,0.8],"linestyle",":","k");
plot([1950,2000],[1.0,1.0],"linestyle",":","k");
plot([1975,2000],[1.2,1.2],"linestyle",":","k");
title("The NYT figure without the splicing");
axis([1400 2000 -1.0 1.3])
set(gca, "yaxislocation", "right","xtick",[1400:50:1950 1997],"ytick",-0.8:0.2:1.2);
subplot(2,1,2); hold on;
plot(MBH98(I98_pre_instru,1),MBH98(I98_pre_instru,2)*1.8,"linestyle","-","linewidth",1,'k');
plot(MBH99_raw(I98_instru,1),MBH99_raw(I98_instru,2)*1.8,"linestyle","-","linewidth",1,"k");
plot([1400,2000],[0,0],"linestyle","-","k");
plot([1400,2000],[0.2,0.2],"linestyle",":","k");
plot([1400,2000],[0.4,0.4],"linestyle",":","k");
plot([1400,2000],[-0.2,-0.2],"linestyle",":","k");
plot([1400,2000],[-0.4,-0.4],"linestyle",":","k");
plot([1400,2000],[-0.6,-0.6],"linestyle",":","k");
plot([1400,2000],[-0.8,-0.8],"linestyle",":","k");
plot([1900,2000],[0.6,0.6],"linestyle",":","k");
plot([1925,2000],[0.8,0.8],"linestyle",":","k");
plot([1950,2000],[1.0,1.0],"linestyle",":","k");
plot([1975,2000],[1.2,1.2],"linestyle",":","k");
title("The NYT figure as published (with the splicing)");
axis([1400 2000 -1.0 1.3]);
set(gca, "yaxislocation", "right","xtick",[1400:50:1950 1997],"ytick",-0.8:0.2:1.2);

figure(1); 
print "trick_MBH98fig5b.jpg"
print "trick_MBH98fig5b.png"
print "trick_MBH98fig5b.pdf"
figure(2); 
print "trick_MBH99fig3a.jpg"
print "trick_MBH99fig3a.png"
print "trick_MBH99fig3a.pdf"
figure(3); 
print "trick_MBH98PRfig2.jpg"
print "trick_MBH98PRfig2.png"
print "trick_MBH98PRfig2.pdf"
figure(4); 
print "trick_MBH98NYTfig.jpg"
print "trick_MBH98NYTfig.png"
print "trick_MBH98NYTfig.pdf"
