Initial commit, added matlab files

2Drebin.m

@@ -0,0 +1,47 @@
+function [binMean,binStandardDev]=2Drebin( x,y,z,sx,sy )
+
+
+mask=~(isnan(x) | isnan(y) | isnan(z));
+x=x(mask);
+y=y(mask);
+z=z(mask)
+
+XbotEdge=min(x(:));
+XtopEdge=max(x(:));
+
+YbotEdge=min(x(:));
+YtopEdge=max(x(:));
+
+XbinEdges=XbotEdge:sx:XtopEdge;
+YbinEdges=YbotEdge:sy:YtopEdge;
+
+[~,XwhichBin]=histc(x,XbinEdges);
+[~,YwhichBin]=histc(y,YbinEdges);
+
+% progressbar(0)
+
+for i=0:numel(XbinEdges)
+    for j=0:numel(YbinEdges)
+
+    XflagBinMembers=(XwhichBin==i);
+    YflagBinMembers=(YwhichBin==i);
+
+    binMembers=z(XflagBinMembers & YflagBinMembers);
+
+    binMean(i+1,j+1)=mean(binMembers);  
+    binStandardDev(i+1,j+1)=std(binMembers);  
+
+%     progressbar(i/numel(binEdges));
+
+    end   
+end
+% 
+% progressbar(1)
+% 
+% binMean=binMean(2:end);
+% binStandardDev=binStandardDev(2:end);
+% binCenters=binEdges+s/2;
+
+
+
+

AddZeroHarm.m

@@ -0,0 +1,6 @@
+function lmcosi=AddZeroHarm(lmcosi,Value)
+
+s=size(lmcosi);
+first_line=zeros(1,s(2));
+first_line(3)=Value;
+lmcosi=[first_line;lmcosi];

AdmittancePCA.m

@@ -0,0 +1,142 @@
+%% PCA of Admittance
+% [coeff,score,latent,tsquared,explained,mu] = pca(Z_t);
+
+ %Z_t=density_t;
+
+ sz=size(Z_t);
+
+[w,pc,ev,tsquared] = princomp(Z_t);
+
+mean_Z_t=mean(Z_t);
+mean_Z_t=repmat(mean_Z_t,[sz(1) 1]);
+
+Ncomp=3;
+
+RGBcomp=[1 2 3];
+MaxDegreeExp=24;
+% Npc=5;
+for Npc=1:Ncomp
+% Y=w(:,Npc);
+Y=w(:,Npc)';
+
+lmcosi_Y=xyz2plm(Y,MaxDegreeExp,'irr',fii,lambdai);
+[Y_sh(:,:,Npc),lonsh,latsh]=plm2xyz(lmcosi_Y,Resolution);
+[lonsh,latsh]=meshgrid(lonsh,latsh);
+
+AGUaxes
+hp=pcolorm(latsh/180*pi,lonsh/180*pi,Y_sh(:,:,Npc)); shading interp
+cb1=colorbar('FontSize',20);
+% ylabel(cb1,'Effective density [kg/m ^3] ','FontSize',20);
+ %PlotVestaFeatures
+title(['Principal component # ' num2str(Npc)],'FontSize',20);
+
+end
+
+s_sh=size(Y_sh);
+RGB=zeros(s_sh(1),s_sh(2),3);
+RGB(:,:,1)=Y_sh(:,:,RGBcomp(1));
+RGB(:,:,2)=Y_sh(:,:,RGBcomp(2));
+RGB(:,:,3)=Y_sh(:,:,RGBcomp(3));
+
+Rrange= max(max(RGB(:,:,1)))-min(min(RGB(:,:,1)));
+Grange= max(max(RGB(:,:,2)))-min(min(RGB(:,:,2)));
+Brange= max(max(RGB(:,:,3)))-min(min(RGB(:,:,3)));
+
+RGB(:,:,1)=(RGB(:,:,1)-min(min(RGB(:,:,1))))/Rrange;
+RGB(:,:,2)=(RGB(:,:,2)-min(min(RGB(:,:,2))))/Grange;
+RGB(:,:,3)=(RGB(:,:,3)-min(min(RGB(:,:,3))))/Brange;
+
+
+ AGUaxes
+% axesm('mollweid','frame','on','FontSize',24);
+pl_pcr=surfm(latsh/180*pi,lonsh/180*pi,ri); shading interp
+% colormap jet
+
+% plot_pc=surflsrm(latsh/180*pi,lonsh/180*pi,ri);
+%  
+set(pl_pcr,'CData',RGB);
+%PlotVestaFeatures
+
+ShowVestaLow
+lighting phong
+
+light_handle1=light('Style','infinite');
+light_handle2=light('Style','infinite');
+
+
+set(h,'CData',RGB);
+lambda_light=0;
+fi_light=90;
+[x_light, y_light, z_light]=sph2cart(lambda_light/180*pi,fi_light/180*pi,1);
+set(light_handle1,'Position',[x_light, y_light, z_light]);
+
+lambda_light=0;
+fi_light=-90;
+[x_light, y_light, z_light]=sph2cart(lambda_light/180*pi,fi_light/180*pi,1);
+set(light_handle2,'Position',[x_light, y_light, z_light]);
+%% plot principal componet
+take_comp=[1 2 3 4 5];
+
+figure; hold on;
+set(gca,'FontSize',20);
+plot(GoodDegrees,pc(:,take_comp));
+% legend(take_comp)
+
+%% Reconstruct data
+take_comp=[1 2 3 4 5 6];
+
+Z_tr=pc(:,take_comp)*w(:,take_comp)';
+Z_tr=Z_tr+mean_Z_t;
+
+for i=1:sz(2)
+%     density_mean_tr(i)=mean((Z_t(:,i))./Z_t_shape(:,i)*rho);
+%     density_mean_tr(i)=mean(Z_tr(:,i)./Z_t_shape(:,i)*rho);
+%     p_tr(i,:)=polyfit(GoodDegrees,(Z_tr(:,i)./Z_t_shape(:,i)*rho)',1);
+     density_mean_tr(i)=mean(Z_tr(:,i));
+     p_tr(i,:)=polyfit(GoodDegrees,(Z_tr(:,i))',1);
+%      density_std(i)=std(Z_tr(:,i));
+
+end
+
+%
+lmcosi_den_tr=xyz2plm(density_mean_tr,MaxDegreeExp,'irr',fii,lambdai);
+[density_mean_sh_tr,lonsh,latsh]=plm2xyz(lmcosi_den_tr,Resolution);
+[lonsh,latsh]=meshgrid(lonsh,latsh);
+
+AGUaxes
+pcolorm(latsh/180*pi,lonsh/180*pi,density_mean_sh_tr); shading interp
+cb1=colorbar('FontSize',20);
+ylabel(cb1,'Effective density [kg/m ^3] ','FontSize',20);
+PlotVestaFeatures
+title('Effective density ','FontSize',20);
+%
+lmcosi_p1_tr=xyz2plm(p_tr(:,1),MaxDegreeExp,'irr',fii,lambdai);
+[p1_sh_tr,lonsh,latsh]=plm2xyz(lmcosi_p1_tr,Resolution);
+[lonsh,latsh]=meshgrid(lonsh,latsh);
+
+AGUaxes
+pcolorm(latsh/180*pi,lonsh/180*pi,p1_sh_tr); shading interp
+cb1=colorbar('FontSize',20);
+ylabel(cb1,'Effective density slope [kg/m ^3/degree] ','FontSize',20);
+PlotVestaFeatures
+title('Effective density slope ','FontSize',20);
+
+%
+% lmcosi_std_tr=xyz2plm(density_std,MaxDegreeExp,'irr',fii,lambdai);
+% [density_std_sh_tr,lonsh,latsh]=plm2xyz(lmcosi_std_tr,Resolution);
+% [lonsh,latsh]=meshgrid(lonsh,latsh);
+% 
+% AGUaxes
+% pcolorm(latsh/180*pi,lonsh/180*pi,density_std_sh_tr); shading interp
+% cb1=colorbar('FontSize',20);
+% ylabel(cb1,'Effective std [kg/m ^3] ','FontSize',20);
+% PlotVestaFeatures
+% title('Effective std ','FontSize',20);
+
+100*ev(take_comp)/sum(ev)
+
+figure
+
+bar(100*ev(take_comp)/sum(ev))
+
+

AdmittanceTest.m

@@ -0,0 +1,144 @@
+ccc
+
+
+%% Admittance
+
+
+GravityFileName='/Users/antonermakov/Dawn/Gravity/VESTA20G/JGV20G02.SHA';
+[lmcosi_grav,Rref,mu,mu_std]=ReadGRAILGravityModel(GravityFileName);
+
+lmcosi_grav=AddZeroHarm(lmcosi_grav,1);
+
+Resolution=1;
+
+MaxTopoPower=10;
+MaxDegreeTopo=100;
+MaxDegreeGrav=20;
+
+% int_str_fig=figure('color',[0 0 0]);
+% plotplm(lmcosi_mantle_shape_new,[],[],2,1,[],[],[]);
+% alpha(0.5)
+% hold on
+%  
+load VestaHASTALAVESTAshape_sh720.mat
+lmcosi_shape=TruncateGravityModel(lmcosi_shape,MaxDegreeTopo,1);
+
+MeanRadius=lmcosi_shape(1,3);
+
+% lmcosi_shape(:,3)=lmcosi_shape(:,3)/MeanRadius;
+% lmcosi_shape(:,4)=lmcosi_shape(:,4)/MeanRadius;
+
+%% Computing admittance
+
+
+% Adm=SphericalHarmonicAdmittance(lmcosi_grav,lmcosi_shape,'r');
+% 
+% hold on;
+% 
+% n=1:MaxDegreeTopo;
+% 
+% CompDepth=10000;
+% 
+% zeta=CompDepth/MeanRadius;
+% 
+% rho_ratio=3000/3457;
+% 
+% AdmTheor=3./(2*n+1).*((1-zeta).^n).*rho_ratio;
+% 
+% 
+%ri=plm2xyz(lmcosi_shape,1)-MeanRadius;
+% 
+%ri_power=ri.^2;
+% 
+%lmcosi_shape_power=xyz2plm(ri_power,MaxDegreeTopo,'im',[],[],[]);  
+%  
+% 
+% 
+% plot(n,AdmTheor);
+% 
+% set(gca,'YScale','Log');
+% set(gca,'XScale','Log');
+% 
+% xlim([2 20]);
+% ylim([0 0.6]);
+
+%% Admittance of a non-spherical body
+
+Rref=265000;
+mu=17.29e9;
+rho=3457;
+
+
+lmcosi_shape(2:end,3)=lmcosi_shape(2:end,3);
+lmcosi_shape(2:end,4)=lmcosi_shape(2:end,4);
+
+[ri,~,~,~]=plm2xyz(lmcosi_shape,Resolution);
+
+lmcosi_shape(:,3)=lmcosi_shape(:,3)/1000;
+lmcosi_shape(:,4)=lmcosi_shape(:,4)/1000;
+
+lmcosi_shape(:,3)=lmcosi_shape(:,3);
+lmcosi_shape(:,4)=lmcosi_shape(:,4);
+
+lmcosi_grav_shape=TopoSH2GravitySH(ri,mu,rho,Rref,MaxDegreeTopo,MaxDegreeGrav,MaxTopoPower);
+
+% lmcosi_calc=AddZeroHarm(lmcosi_calc,1);
+
+% 
+lmcosi_grav_shape(:,3)=lmcosi_grav_shape(:,3).*(lmcosi_grav_shape(:,1)+1)*mu/(Rref^2)*1e5;
+lmcosi_grav_shape(:,4)=lmcosi_grav_shape(:,4).*(lmcosi_grav_shape(:,1)+1)*mu/(Rref^2)*1e5;
+
+lmcosi_grav(:,3)=lmcosi_grav(:,3).*(lmcosi_grav(:,1)+1)*mu/(Rref^2)*1e5;
+lmcosi_grav(:,4)=lmcosi_grav(:,4).*(lmcosi_grav(:,1)+1)*mu/(Rref^2)*1e5;
+% 
+% lmcosi_calc(:,3)=lmcosi_calc(:,3).*(lmcosi_calc(:,1)+1)*mu/(Rref^2)*1e5;
+% lmcosi_calc(:,4)=lmcosi_calc(:,4).*(lmcosi_calc(:,1)+1)*mu/(Rref^2)*1e5;
+
+Adm_homo=SphericalHarmonicAdmittance(lmcosi_grav_shape,lmcosi_shape);
+Adm_obs=SphericalHarmonicAdmittance(lmcosi_grav,lmcosi_shape);
+
+
+% Adm_calc=SphericalHarmonicAdmittance(lmcosi_calc,lmcosi_shape);
+
+
+% set(gca,'YScale','Log');
+% set(gca,'XScale','Log');
+
+
+n=1:MaxDegreeGrav;
+
+hold on;
+
+plot(n,Adm_homo,'o-','Color','r','LineWidth',4,'MarkerSize',2);
+plot(n,Adm_obs,'o-','Color','b','LineWidth',4,'MarkerSize',2);
+
+
+
+title('Gravity/Topography Admittance','FontSize',25,'FontName','Times');
+xlabel('Degree','FontSize',25,'FontName','Times');
+
+grid on;
+hold on;
+
+AdmTheor=3./(2*(n)+1).*((MeanRadius/Rref).^n).*(n+1)*mu/(Rref^2)*1e5/1000;
+
+plot(n,AdmTheor,'Color','k','LineWidth',4,'MarkerSize',2);
+
+
+
+% plot(1:16,Adm_calc,'o-','Color','g','LineWidth',4,'MarkerSize',2);
+
+
+% ylim([0 1]);
+
+%%
+
+legend({'homogenous','observed','linear','model'},'fontsize',20);
+
+
+set(gca,'FontSize',20)
+
+% set(gca,'YScale','Log');
+
+% ylim([0.04,1]);
+

AiryAdmittance.m

@@ -0,0 +1,9 @@
+function Z=AiryAdmittance(R,g,nu,rhomean,rhomantle,rhocrust,E,d,n,tl)
+
+l=CompensationLevel(R,g,nu,rhomean,rhomantle,rhocrust,E,d,n);
+
+% l=1;
+
+Coef=1.6*1e5/R*1000;
+
+Z=Coef*3./(2.*n+1).*(rhocrust./rhomean).*(1-((1-tl./R).^n).*l).*(n+1);

AiryAdmittanceFun.m

@@ -0,0 +1,14 @@
+function Z=AiryAdmittanceFun(x,n)
+
+g=1.6;
+nu=0.25;
+
+% rhocrust=2600;
+rhomantle=3220;
+rhomean=3344;
+Rref=1738000;
+
+% h_comp=2*x(2);
+% (R,g,nu,rhomean,rhomantle,rhocrust,E,d,n,tl)
+Z=AiryAdmittance(Rref,g,nu,rhomean,rhomantle,x(3),x(1),x(2),n,2*x(2));
+

AxesOrientation.m

@@ -0,0 +1,16 @@
+function [lambda,fi]=AxesOrientation(evecs)
+
+a_vec=evecs(:,1);
+b_vec=evecs(:,2);
+c_vec=evecs(:,3);
+
+lambda_a=atan2(a_vec(2),a_vec(1))*180/pi;
+lambda_b=atan2(b_vec(2),b_vec(1))*180/pi;
+lambda_c=atan2(c_vec(2),c_vec(1))*180/pi;
+
+fi_a=atan2(a_vec(3),norm([a_vec(1) a_vec(2)]))*180/pi;
+fi_b=atan2(b_vec(3),norm([b_vec(1) b_vec(2)]))*180/pi;
+fi_c=atan2(c_vec(3),norm([c_vec(1) c_vec(2)]))*180/pi;
+
+lambda=[lambda_a lambda_b lambda_c];
+fi=[fi_a fi_b fi_c];