% path(path,'/globec/mcgillic/matlab')
path(path,'/globec/mcgillic/matlab_cen')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FORWARD AND INVERSE SOLUTION VELOCITY PLOTS
%
% Load .m3d, .dat, .s2c and grid data
%
runtype=input('Enter the runtype (fsvd, svd or wls): ','s')
%
file_input_m3d='FUNDY5.M2.m3d'
file_output_dat=['M2_',runtype,'.dat']
file_output_s2c=['M2_',runtype,'_zeta.s2c']
file_output_.bcs.s2c=['M2_',runtype,'_zeta.bcs.s2c']
[OBS_M3D,gridname,year,ncol,gmt]=read_obs(file_input_m3d);
[INV_M3D]=loadname(file_output_dat);
[INV_S2C,freq,gridname]=read_s2c(file_output_s2c);
amp=INV_S2C(:,2);
pha=INV_S2C(:,3);
[INV_S2C,freq,gridname]=read_s2c(file_output_.bcs.s2c);
[real,aimag]=ap_to_ri_deg(INV_S2C(:,2),INV_S2C(:,3));
%
gridname=blank(gridname(1:length(gridname)-1));
gridstruct=loadgrid(gridname);
in=gridstruct.e;x=gridstruct.x;y=gridstruct.y;z=gridstruct.z;
bnd=gridstruct.bnd;
%
% Report statistics
%
V=OBS_M3D;
rms(V(:,4)+sqrt(-1.0)*V(:,5))
V=INV_M3D;
rms(V(:,4)+sqrt(-1.0)*V(:,5))
rms(V(:,6)+sqrt(-1.0)*V(:,7))
%
% Plot results
%
landscapefigure
whitebg('w')
%
%
subplot(2,2,1)
plotbnd(gridstruct)
hold on
V=OBS_M3D;
hv=vecplot2(V(:,2),V(:,3),V(:,4),V(:,5),1.0,'m/s',175000.,-25000.);
set(hv,'color','blue')
title(verbatim([file_input_m3d,' Observations: RMS = ',num2str(rms(V(:,4)+sqrt(-1.0)*V(:,5))),' m/s']))
axis('tight')
%
%
subplot(2,2,2)
plotbnd(gridstruct)
hold on
V=INV_M3D;
hv=vecplot2(V(:,2),V(:,3),V(:,4),V(:,5),1.0,'m/s',175000.,-25000.);
set(hv,'color','red')
title(verbatim([file_output_dat,' Inversion: RMS = ',num2str(rms(V(:,4)+sqrt(-1.0)*V(:,5))),' m/s']))
axis('tight')
%
%
subplot(2,2,3)
plotbnd(gridstruct)
hold on
V=INV_M3D;
hv=vecplot2(V(:,2),V(:,3),V(:,6),V(:,7),0.1,'m/s',175000.,-25000.);
set(hv,'color','g')
axis('tight')
%
scalar=amp;
smin=0.0;smax=1.4;ibw=0.2;
[smin,smax,ibw]=colorband_fe(in,x,y,bnd,scalar,smin,smax,ibw)
title(['Inverse Error: RMS = ',num2str(rms(V(:,6)+sqrt(-1.0)*V(:,7))),' m/s'])
xlabel('Tidal Amplitude (m)')
axis('tight')
%
%
subplot(2,2,4)
plotbnd(gridstruct)
hold on
V=INV_M3D;
hv=vecplot2(V(:,2),V(:,3),V(:,6),V(:,7),0.1,'m/s',175000.,-25000.);
set(hv,'color','g')
%
scalar=pha;
smin=120.0;smax=330.0;ibw=30.0;
[smin,smax,ibw]=colorband_fe(in,x,y,bnd,scalar,smin,smax,ibw)
title(['Inverse Error: RMS = ',num2str(rms(V(:,6)+sqrt(-1.0)*V(:,7))),' m/s'])
xlabel('Tidal Phaselag (deg gmt)')
axis('tight')
%
%
%
figure
hold on
plot(real,'b-');plot(aimag,'r--');
title('Real and imaginary parts of BCS')
ylabel('elevation (m)')
xlabel('boundary node number')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
