% Matlab routine to load netcdf sst files and create a VRML
% Elevation Grid for inclusion in the TAO demo world.
% Matlab must have the netcdf tools installed, 
% see http://crusty.er.usgs.gov/mexcdf.html if needed
% Created 9/28/98 by C. Moore
clear
for n=1:11
n
%
% open netcdf file - note: the netcdf tools must be available to
% Matlab, and ncstartup.m must be run before using this script
%
indir='../dyn/lat_lon/';
outdir='../../public_html/';
if n<10
 filei=[indir,'980',int2str(n),'_dynht_xy_lf.cdf'];
 fileo=[outdir,'980',int2str(n),'_dyn.wrl'];
else
 filei=[indir,'98',int2str(n),'_dynht_xy_lf.cdf'];
 fileo=[outdir,'98',int2str(n),'_dyn.wrl'];
end
f=netcdf(filei,'nowrite');
data=f{6}(:);
%sst=tsst(1,1,:,:);
lat=f{4}(:);
lon=f{5}(:);
ncclose;
%
% Open the output file, and print the header
%
fid=fopen(fileo,'w');
fprintf(fid,'#VRML V2.0 utf8\n');
fprintf(fid,'# created by mkdyn.m Matlab routine\n');
fprintf(fid,'# Chris Moore UW/NOAA\n');
fprintf(fid,'Transform {\n');
% Watch out!  My world might need a translation or rotation!
fprintf(fid,'rotation 1 0 0 3.14159\n');
fprintf(fid,'children [\n');
fprintf(fid,'Shape {\n');
fprintf(fid,'  appearance Appearance {\n');
fprintf(fid,'    texture ImageTexture {url "sst%d.jpeg"}\n',n);
fprintf(fid,'    textureTransform TextureTransform {\n');
fprintf(fid,'      center .5 .5\n');
fprintf(fid,'      translation .315 5.8\n');
fprintf(fid,'      scale 0.775 .15}\n');
fprintf(fid,'  }\n');
fprintf(fid,'  geometry ElevationGrid {\n');
fprintf(fid,'	solid FALSE\n');
% The next five lines can be generalized to arbitrary Dimension
fprintf(fid,'     xDimension 129\n');
fprintf(fid,'     zDimension 19\n');
fprintf(fid,'     xSpacing 4.0\n');
fprintf(fid,'     zSpacing 4.0\n');
fprintf(fid,'     creaseAngle 0.785\n');
fprintf(fid,'     height [\n');
%
% Demean and scale the data to create height
%
SCALE = 1.0;
OFFSET = -20.0;
BADDATA = 1.0e5;
i=find(data<BADDATA);
mn=sum(data(i))/length(i);
datadm = -1*(data-mn);
height = SCALE.*datadm+OFFSET;
for i=1:19
 for j=1:129
  if data(i,j)<BADDATA
   fprintf(fid,' %1.2f,', height(i,j));
  else
   fprintf(fid,' 0.00,');
  end
 end
 fprintf(fid,'\n');
end
%
% Print more text header
%
fprintf(fid,'     ]\n');
fprintf(fid,'     colorPerVertex TRUE\n');
fprintf(fid,'     color Color {\n');
fprintf(fid,'       color [\n');
%
% the data are scaled to fit our colortable
%
load colortable
mincolor=1;
maxcolor=length(colortable)-1;
maxdata=150;
mindata=80;
m=(maxcolor-mincolor)/(maxdata-mindata);
b=1-mindata*m;
for i=1:19,
 for j=1:129,
  ind=floor(m*data(i,j)+b);
   if (ind > mincolor)&(ind < maxcolor+1)
    fprintf(fid,' %1.4f %1.4f %1.4f,',colortable(ind,:));
   else
    fprintf(fid,' %1.4f %1.4f %1.4f,',colortable(length(colortable),:));
   end
 end
 fprintf(fid,'\n');
end
%
% Print the closing text
%
fprintf(fid,'       ]\n');
fprintf(fid,'     }\n');
fprintf(fid,'  }\n');
fprintf(fid,'}\n');
fprintf(fid,']\n');
fprintf(fid,'}\n');
fclose(fid);
%
% uncomment the following to contour each to view:
%
%v=20:30; % use mindata:maxdata?
%[cs,h]=contourf(lon,lat,data,v);
%clabel(cs,h,'manual');
%axis('image');
end