Index: ../trunk-jpl/examples/EsaGRACE/grace.m =================================================================== --- ../trunk-jpl/examples/EsaGRACE/grace.m (revision 22809) +++ ../trunk-jpl/examples/EsaGRACE/grace.m (nonexistent) @@ -1,121 +0,0 @@ -% map grace loads in meters [m] of water equivalent height -function water_load = grace(index,lat,lon,tmin,tmax); - - %compute centroids using (lat,lon) data {{{ - ne = length(index); % number of elements - % lat -> [0,180]; long -> [0,360] to compute centroids - lat=90-lat; - lon(lon<0)=180+(180+lon(lon<0)); - - ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); - bx_0=lat(index(:,2)); by_0=lon(index(:,2)); - cx_0=lat(index(:,3)); cy_0=lon(index(:,3)); - % find whether long is 0 or 360! This is important to compute centroids as well as elemental area - for ii=1:ne - 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; - % - lat_element=(ax_0+bx_0+cx_0)/3; - lon_element=(ay_0+by_0+cy_0)/3; - - % [lat,lon] \in [-90:90,-180,180]; - lat_element_0 = 90-lat_element; lon_element_0 = lon_element; - lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180; - % }}} - - % Monthly GRACE data - filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; - time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC - long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75 - lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75 - rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm] - - time_yr = 2002+time_0/365; % [yr] - - [nn_0,mm_0] = size(squeeze(rec(:,:,1))); - for jj=1:mm_0 % chose a latitude - for kk=1:nn_0 - ii=nn_0*(jj-1)+kk; - lat_0(ii)=lati_0(jj); - lon_0(ii)=long_0(kk); - tws_monthly(:,ii) = rec(kk,jj,:); - end - end - % - %% translate variables as grace-related variables -- so I do not need to do too much editing - grace_monthly=tws_monthly; - grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0; - - % detrend over the entire time period - grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column! - - % fill out the blanks {{{ - - lat_grace=lat_0; - lon_grace=lon_0; - num_org=length(lon_grace); - - qq=1; mm=1; - for jj=2:num_org-1 - if (lat_grace(jj)~=lat_grace(jj+1)) - lat_new(qq)=lat_grace(jj); - lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1)); - load_new(:,qq)=grace_monthly(:,mm); - lat_new(qq+1)=lat_grace(jj); - lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1)); - load_new(:,qq+1)=grace_monthly(:,jj); - qq=qq+2; - mm=jj+1; % to find out the value for monthly data - end - end - - num_add=length(lat_new); - num_plus=num_org+num_add; - - lat_grace_plus=zeros(num_plus,1); - lon_grace_plus=zeros(num_plus,1); - load_grace_plus=zeros(length(time_0),num_plus); - - lat_grace_plus(1:num_org)=lat_grace; - lat_grace_plus(1+num_org:num_plus)=lat_new; - lon_grace_plus(1:num_org)=lon_grace; - lon_grace_plus(1+num_org:num_plus)=lon_new; - load_grace_plus(:,1:num_org)=grace_monthly; - load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:); - % }}} - - % choose selected months ONLY - [diff1,pos1] = min(abs(tmin-time_yr)); - [diff2,pos2] = min(abs(tmax-time_yr)); - - time_yr=time_yr(pos1:pos2); - load_grace_plus=load_grace_plus(pos1:pos2,:); - num_yr=length(time_yr); - water_load_0=zeros(ne,num_yr); - - for jj=1:num_yr - water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element); - disp([num2str(jj),' of ',num2str(num_yr),' months done!']); - end - - water_load=water_load_0/100; % cm -> meters of water - water_load(isnan(water_load))=0; - Property changes on: ../trunk-jpl/examples/EsaGRACE/grace.m ___________________________________________________________________ Deleted: svn:executable ## -1 +0,0 ## -* \ No newline at end of property Index: ../trunk-jpl/examples/EsaWahr/wahr.m =================================================================== --- ../trunk-jpl/examples/EsaWahr/wahr.m (revision 22809) +++ ../trunk-jpl/examples/EsaWahr/wahr.m (nonexistent) @@ -1,46 +0,0 @@ -% compute semi-analytic solutions for a disc load -function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); - - disc_rad = disc_rad/1000; % km - % compute P(x), dP(x)/dx, d2P(x)/dx2 - %--------------------------------------------------------------------- - % compute p_value - theta=km2deg(xi/1000)'; - ang = theta/180*pi; - alpha=cos(ang); - m=length(alpha); - n=length(love_h)-1; - p_value = p_polynomial_value(m,n,alpha); - p_prime = p_polynomial_prime(m,n,alpha); - %--------------------------------------------------------------------- - nn=[0:n]; - nn_plus_1=nn+1; - - % disc radius in degree - disc_rad = km2deg(disc_rad)/180*pi; - tau=zeros(size(love_h)); - tau(1) = 0.5*(1-cos(disc_rad)); % tau_0 - p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad)); - p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad)); - for jj=2:n+1 - nnn = jj-1; - tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); - end - - const=zeros(size(love_h)); - for jj=1:n+1 - const(jj) = 1/(2*(jj-1)+1); - end - - disc=sum(bsxfun(@times,p_value,tau'),2); - - g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2); - g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2); - - % coeff - coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; - - % vertical and horizontal solutions in mm - vert = g1*coeff*1000; % mm - horz = g5*coeff*1000; % mm - Property changes on: ../trunk-jpl/examples/EsaWahr/wahr.m ___________________________________________________________________ Deleted: svn:executable ## -1 +0,0 ## -* \ No newline at end of property