source: issm/oecreview/Archive/24684-25833/ISSM-25401-25402.diff

../trunk-jpl/src/m/contrib/adhikari/latelonge.m

@@ +1 @@
+function [late,longe] = latelonge(md_mesh) 
+%latelonge :: computes lat,long at elemental centroids on spherical surface
+%
+%   Usage:
+%      [late,longe]=latelonge(md.mesh)
+
+nv = md_mesh.numberofvertices; 
+if length(md_mesh.lat)~=nv | length(md_mesh.long)~=nv 
+        error('lat,long not defined properly at vertices.');
+end
+
+% lat -> [0,180]; long -> [0,360] to compute centroids 
+lat=90-md_mesh.lat;             lon=md_mesh.long; 
+lon(lon<0)=180+(180+lon(lon<0)); 
+
+ax_0=lat(md_mesh.elements(:,1)); ay_0=lon(md_mesh.elements(:,1)); 
+bx_0=lat(md_mesh.elements(:,2)); by_0=lon(md_mesh.elements(:,2)); 
+cx_0=lat(md_mesh.elements(:,3)); cy_0=lon(md_mesh.elements(:,3)); 
+% find whether long is 0 or 360! This is important to compute centroids as well as elemental area 
+for ii=1:md_mesh.numberofelements
+        if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
+                if ay_0(ii)==0
+                        ay_0(ii)=360;
+                end 
+                if by_0(ii)==0
+                        by_0(ii)=360; 
+                end 
+                if cy_0(ii)==0 
+                        cy_0(ii)=360; 
+                end
+        end 
+end
+% correction at the north pole 
+ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; 
+by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; 
+cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; 
+% correction at the south pole 
+ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; 
+by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; 
+cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; 
+% 
+late=(ax_0+bx_0+cx_0)/3; 
+longe=(ay_0+by_0+cy_0)/3;
+
+% bak to [-90 90] [-180 180] ranges. 
+late = 90-late; 
+longe(longe>180) = longe(longe>180)-360; 
+
+