% PS05_SEMIFIT an m-file to fit theoretical semivariograms to supplied 
%   empirical semivariances. This program will fit spherical models to 
%   the empirical values. The results are displayed to the screen and 
%   written out to a separate ASCII file. If you want to see the results
%   of using the exponential model, replace SPHERMOD with EXPMOD in the
%   calls to NLLEASQR.
% 
%   NOTE: the results of the exponential model fit are stored in 
%   PS05_SEMIFIT_EXP.ASC and PS05_SEMIFIT_EXP.EPS
%
%   NOTE: each vector fit has its own "stop" length. This is to ensure
%   that no periodic structures are included in the data fed to the 
%   theoretical semivariogram. 
%
% Started 2002:10:21 D. Glover, WHOI
% Modif'd 2008:10:27 DMG changed the name from PS04 to PS05

% Start a stopwatch

t_init_semifit=clock;

% Read in the supplied empirical semivariances

load semivar;

% Compute lag vector in degrees for comparison to images

lagdeg=[1:113]*0.088496;	% 0.088496 deg/pixel

% Set pin and fit spherical model

pin=[0.1,0.5,1];

nstop11x=24;		% length of fitted semivariogram vector

[f11x,p11x,kvg,iter,corp,covp,covr,stdresid,Z,r2x11] = nlleasqr(lagdeg(1:nstop11x)',S11x(1:nstop11x)',pin,'sphermod');
p11xerr=sqrt(diag(covp));

nstop11y=25;		% length of fitted semivariogram vector

[f11y,p11y,kvg,iter,corp,covp,covr,stdresid,Z,r2y11] = nlleasqr(lagdeg(1:nstop11y)',S11y(1:nstop11y)',pin,'sphermod');
p11yerr=sqrt(diag(covp));

nstop12x=24;		% length of fitted semivariogram vector

[f12x,p12x,kvg,iter,corp,covp,covr,stdresid,Z,r2x21] = nlleasqr(lagdeg(1:nstop12x)',S12x(1:nstop12x)',pin,'sphermod');
p12xerr=sqrt(diag(covp));

nstop12y=34;		% length of fitted semivariogram vector

[f12y,p12y,kvg,iter,corp,covp,covr,stdresid,Z,r2y21] = nlleasqr(lagdeg(1:nstop12y)',S12y(1:nstop12y)',pin,'sphermod');
p12yerr=sqrt(diag(covp));

nstop13x=26;		% length of fitted semivariogram vector

[f13x,p13x,kvg,iter,corp,covp,covr,stdresid,Z,r2x13] = nlleasqr(lagdeg(1:nstop13x)',S13x(1:nstop13x)',pin,'sphermod');
p13xerr=sqrt(diag(covp));

nstop13y=34;		% length of fitted semivariogram vector

[f13y,p13y,kvg,iter,corp,covp,covr,stdresid,Z,r2y13] = nlleasqr(lagdeg(1:nstop13y)',S13y(1:nstop13y)',pin,'sphermod');
p13yerr=sqrt(diag(covp));

nstop22x=56;		% length of fitted semivariogram vector

[f22x,p22x,kvg,iter,corp,covp,covr,stdresid,Z,r2x22] = nlleasqr(lagdeg(1:nstop22x)',S22x(1:nstop22x)',pin,'sphermod');
p22xerr=sqrt(diag(covp));

nstop22y=34;		% length of fitted semivariogram vector

[f22y,p22y,kvg,iter,corp,covp,covr,stdresid,Z,r2y22] = nlleasqr(lagdeg(1:nstop22y)',S22y(1:nstop22y)',pin,'sphermod');
p22yerr=sqrt(diag(covp));

nstop23x=42;		% length of fitted semivariogram vector

[f23x,p23x,kvg,iter,corp,covp,covr,stdresid,Z,r2x23] = nlleasqr(lagdeg(1:nstop23x)',S23x(1:nstop23x)',pin,'sphermod');
p23xerr=sqrt(diag(covp));

nstop23y=46;		% length of fitted semivariogram vector

[f23y,p23y,kvg,iter,corp,covp,covr,stdresid,Z,r2y23] = nlleasqr(lagdeg(1:nstop23y)',S23y(1:nstop23y)',pin,'sphermod');
p23yerr=sqrt(diag(covp));

nstop33x=28;		% length of fitted semivariogram vector

[f33x,p33x,kvg,iter,corp,covp,covr,stdresid,Z,r2x33] = nlleasqr(lagdeg(1:nstop33x)',S33x(1:nstop33x)',pin,'sphermod');
p33xerr=sqrt(diag(covp));

nstop33y=32;		% length of fitted semivariogram vector

[f33y,p33y,kvg,iter,corp,covp,covr,stdresid,Z,r2y33] = nlleasqr(lagdeg(1:nstop33y)',S33y(1:nstop33y)',pin,'sphermod');
p33yerr=sqrt(diag(covp));

% Plot up the raw N-S/E-W semivariances and their fitted semivariograms

figure

subplot(321);plot(lagdeg,S11x,'r+');hold
plot(lagdeg,S11y,'bo')
plot(lagdeg(1:nstop11x),f11x(1:nstop11x),'r','linewidth',2)
plot(lagdeg(1:nstop11y),f11y(1:nstop11y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit',4)
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 130 vs. 130')

subplot(322);plot(lagdeg,S12x,'r+');hold
plot(lagdeg,S12y,'bo')
plot(lagdeg(1:nstop12x),f12x(1:nstop12x),'r','linewidth',2)
plot(lagdeg(1:nstop12y),f12y(1:nstop12y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit')
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 130 vs. 131')

subplot(323);plot(lagdeg,S22x,'r+');hold
plot(lagdeg,S22y,'bo')
plot(lagdeg(1:nstop22x),f22x(1:nstop22x),'r','linewidth',2)
plot(lagdeg(1:nstop22y),f22y(1:nstop22y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit')
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 131 vs. 131')

subplot(324);plot(lagdeg,S13x,'r+');hold
plot(lagdeg,S13y,'bo')
plot(lagdeg(1:nstop13x),f13x(1:nstop13x),'r','linewidth',2)
plot(lagdeg(1:nstop13y),f13y(1:nstop13y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit')
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 130 vs. 132')

subplot(325);plot(lagdeg,S33x,'r+');hold
plot(lagdeg,S33y,'bo')
plot(lagdeg(1:nstop33x),f33x(1:nstop33x),'r','linewidth',2)
plot(lagdeg(1:nstop33y),f33y(1:nstop33y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit')
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 132 vs. 132')

subplot(326);plot(lagdeg,S23x,'r+');hold
plot(lagdeg,S23y,'bo')
plot(lagdeg(1:nstop23x),f23x(1:nstop23x),'r','linewidth',2)
plot(lagdeg(1:nstop23y),f23y(1:nstop23y),'b','linewidth',2)
grid;
axis([0 5 0 0.6001]);
legend('E-W empirical','N-S empirical','E-W fit','N-S fit')
ylabel('Semivariance (^{\circ}C^{2})')
xlabel('Lag (degrees)')
title('AVHRR 1987 day 131 vs. 132')
text(4,-0.1,'PS05c 2002');

% Display the fitted parameters and their errors

format compact  % Put into compact format to reduce number of carriage returns

disp('Parameters (nugget, sill, range)');
disp(' ');
disp('     P11x      P11y      P22x      P22y      P33x      P33y');
[p11x p11y p22x p22y p33x p33y]
disp(' ');
disp('     P12x      P12y      P13x      P13y      P23x      P23y');
[p12x p12y p13x p13y p23x p23y]
disp(' ');

disp('Parameter errors (nugget, sill, range)')
disp(' ');
disp('   err11x    err11y    err22x    err22y    err33x    err33y');
[p11xerr p11yerr p22xerr p22yerr p33xerr p33yerr]
disp(' ');
disp('   err12x    err12y    err13x    err13y    err23x    err23y');
[p12xerr p12yerr p13xerr p13yerr p23xerr p23yerr]
disp(' ');

% Write results to an ASCII file

semifit_fid=fopen('ps05_semifit.asc','wt+');

count=fprintf(semifit_fid,'%s\n\n','Parameters (nugget, sill, range)');
count=fprintf(semifit_fid,'%s\n','     P11x       P11y       P22x       P22y       P33x       P33y');
count=fprintf(semifit_fid,'%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',[p11x'; p11y'; p22x'; p22y'; p33x'; p33y']);
count=fprintf(semifit_fid,'\n%s\n','     P12x       P12y       P13x       P13y       P23x       P23y');
count=fprintf(semifit_fid,'%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',[p12x'; p12y'; p13x'; p13y'; p23x'; p23y']);
count=fprintf(semifit_fid,'\n%s\n\n','Parameter errors (nugget, sill, range)');
count=fprintf(semifit_fid,'%s\n','   err11x     err11y     err22x     err22y     err33x     err33y');
count=fprintf(semifit_fid,'%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',[p11xerr'; p11yerr'; p22xerr'; p22yerr'; p33xerr'; p33yerr']);
count=fprintf(semifit_fid,'\n%s\n','   err12x     err12y     err13x     err13y     err23x     err23y');
count=fprintf(semifit_fid,'%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',[p12xerr'; p12yerr'; p13xerr'; p13yerr'; p23xerr'; p23yerr']);

fclose('all');  % remember to close the file, or you won't see the
                % latest write
format loose

% Report semivariance execution time

sprintf('%s%g\n','SemiFit Elasped Time: ',etime(clock,t_init_semifit))