Warning, /utils/matlab/cs_grid/latloncap/llc_vec.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6a75a3b8ef Chri*0001 function h=llc_vec(t,xfld,yfld,fspec);
0002 %function h=llc_vecl(t,xfld,yfld,fspec)
0003 % t - rdnctiles data structure
0004 %
0005 % xfld - field to use for x-component
0006 % yfld - field to use for y-component
0007 % fspec - things to control plot region etc...
0008 % fspec.proj is either 'cart' or 'cyl'. 'sphere' with quiver3 is
0009 % too slow to be useful on my laptop.
0010 % fspec.lathi high latitude value for each region in the plot
0011 % fspec.latlo low latitude value for each region in the plot
0012 % fspec.lonhi high longitude value for each region in the plot
0013 % fspec.lonlo low longitude value for each region in the plot
0014 % fspec.az azimuthal view angle
0015 % fspec.el elevation view angle
0016 % fspec.alx scale arrow x-location
0017 % fspec.aly scale arrow y-location
0018 % fspec.qscal quiver scale factor
0019 % fspec.skip arrow skip 1, every arrow, 2 every two arrows etc...
0020 % fspec.almag scale arrow length
0021 % fspec.xlab xlabel
0022 % fspec.ylab ylabel
0023 % Code that does the scale arrow is a bit of a hack - to say the
0024 % least :-).
0025
0026 % Make a vector plot using quiver
0027 clf; h=[];
0028
0029 xvec=[]; yvec=[]; xloc=[]; yloc=[];
0030 xp=[];yp=[];
0031 for tn=1:length(t)
0032 uvec=getfield(t(tn).var,xfld );
0033 vvec=getfield(t(tn).var,yfld );
0034 xc =getfield(t(tn).var,'XC');
0035 yc =getfield(t(tn).var,'YC');
0036 xg =getfield(t(tn).var,'XG');
0037 yg =getfield(t(tn).var,'YG');
0038 % Make list of water and land points
0039 ikeep=[]; iland=[];
0040 i2land=[];j2land=[];
0041 for i=1:length(fspec.lathi)
0042 lathi=fspec.lathi(i);
0043 latlo=fspec.latlo(i);
0044 lonhi=fspec.lonhi(i);
0045 lonlo=fspec.lonlo(i);
0046 ikeep=[ikeep' find( xc >= lonlo ...
0047 & xc <= lonhi ...
0048 & yc >= latlo ...
0049 & yc <= lathi ...
0050 & (uvec.^2+vvec.^2) ~= 0 )']';
0051 iland=[iland' find( xc >= lonlo ...
0052 & xc <= lonhi ...
0053 & yc >= latlo ...
0054 & yc <= lathi ...
0055 & (uvec.^2+vvec.^2) == 0 )']';
0056 end
0057 ikeep=unique(ikeep); iland=unique(iland);
0058 % Extract vector values for water points
0059 xvec=[xvec' uvec(ikeep)']';
0060 yvec=[yvec' vvec(ikeep)']';
0061 xloc=[xloc' xc(ikeep)']';
0062 yloc=[yloc' yc(ikeep)']';
0063 % Create patch array for land points
0064 [tmpi,tmpj]=ind2sub(size(xc),iland);
0065 t(tn).xp=ones(length(tmpi),5);
0066 t(tn).yp=ones(length(tmpi),5);
0067 for l=1:length(tmpi)
0068 t(tn).xp(l,:)=[xg(tmpi(l) , tmpj(l) ) ...
0069 xg(tmpi(l)+1, tmpj(l) ) ...
0070 xg(tmpi(l)+1, tmpj(l)+1) ...
0071 xg(tmpi(l) , tmpj(l)+1) ...
0072 xg(tmpi(l) , tmpj(l) )];
0073 t(tn).yp(l,:)=[yg(tmpi(l) , tmpj(l) ) ...
0074 yg(tmpi(l)+1, tmpj(l) ) ...
0075 yg(tmpi(l)+1, tmpj(l)+1) ...
0076 yg(tmpi(l) , tmpj(l)+1) ...
0077 yg(tmpi(l) , tmpj(l) )];
0078 end
0079 end
0080
0081 if strcmp(fspec.proj,'cyl')
0082 % Do cylindrical plot (would be nice to use quiver3, but its very slow)
0083 [xloc_sp, yloc_sp, zloc_sp]=sph2cart(deg2rad(xloc+360),deg2rad(yloc),6.4e3);
0084 ca=cos(deg2rad(xloc+360)); sa=sin(deg2rad(xloc+360));
0085 if fspec.el == -90
0086 yvec=-yvec;
0087 end
0088 xvec_sp=-sa.*xvec-ca.*yvec; yvec_sp=+ca.*xvec-sa.*yvec;
0089 % Add reference arrow
0090 xloc_sp(end+1)=fspec.alx;
0091 yloc_sp(end+1)=fspec.aly;
0092 xvec_sp(end+1)=fspec.axmag;
0093 yvec_sp(end+1)=fspec.aymag;
0094 hq=quiver(xloc_sp,yloc_sp,xvec_sp,yvec_sp,fspec.scal);
0095 h=[h;hq];
0096 axis equal;axis tight
0097 hold on
0098 for tn=1:length(t)
0099 xp=t(tn).xp; yp=t(tn).yp;
0100 [xp_sp,yp_sp,zp_sp]=sph2cart(deg2rad(xp+360),deg2rad(yp+360),6.4e3);
0101 hp=patch(xp_sp',yp_sp',3);h=[h;hp];
0102 c=[0.6 0.6 0.6];
0103 e=[ 0 0 0];
0104 set(hp,'FaceColor',c,'EdgeColor',e)
0105 set(hp,'EraseMode','background')
0106 end
0107 xlabel('km');
0108 ylabel('km');
0109 view(fspec.az,fspec.el);
0110 % Add label to reference arrow
0111 alab=sprintf('%.1f m s^{-1}',sqrt(fspec.axmag.^2+fspec.aymag.^2));
0112 if fspec.axmag == 0
0113 text(fspec.alx+150, fspec.aly, alab);
0114 else
0115 text(fspec.alx, fspec.aly+150, alab);
0116 end
0117 end
0118
0119 if strcmp(fspec.proj,'cart')
0120 if fspec.az == 180
0121 xloc(find(xloc<0))=xloc(find(xloc<0))+360;
0122 end
0123 % Add reference arrow
0124 xloc(end+1)=fspec.alx;
0125 yloc(end+1)=fspec.aly;
0126 xvec(end+1)=fspec.axmag;
0127 yvec(end+1)=fspec.aymag;
0128 hold on
0129 for tn=1:length(t)
0130 xp=t(tn).xp; yp=t(tn).yp;
0131 if fspec.az == 180
0132 xp(find(xp<0))=xp(find(xp<0))+360;
0133 end
0134 hp=patch(xp',yp',3);h=[h;hp];
0135 c=[0.6 0.6 0.6];
0136 e=[ 0 0 0];
0137 set(hp,'FaceColor',c,'EdgeColor',c)
0138 end
0139 sk=fspec.skip;
0140 hq=quiver(xloc(1:sk:end),yloc(1:sk:end),xvec(1:sk:end),yvec(1:sk:end),fspec.scal);
0141 axis equal;axis tight
0142 h=[h;hq];
0143 xlabel('degrees');
0144 ylabel('degrees');
0145 % Add label to reference arrow
0146 alab=sprintf('%.1f m s^{-1}',sqrt(fspec.axmag.^2+fspec.aymag.^2));
0147 if fspec.axmag == 0
0148 text(fspec.alx+1, fspec.aly, alab);
0149 else
0150 text(fspec.alx, fspec.aly+1, alab);
0151 end
0152 end
0153
0154 return
0155 end
