% Momenten-Tensor
format compact
clear
str=input('Strike? ')/180*pi;  % bei Ungerer S
dip=input('Dip? ')/180*pi;     % bei Ungerer F
rak=input('Rake? ')/180*pi;    % bei Ungerer D

%          M0            Staerke des Momententensors in [dyn*cm]       
%          Mxx,Myy,Mzz   normierter Momententensor  M=(ni*fj+nj*fi)    
%          Mxy,Mxz,Myz   (Bem:seismischer Momententensor MS=M0*M)      
%          Bem2: Die Umrechnung von Streichen(S),Fallen(F) und Dislo-  
%          kationswinkel(D) beim Double-Couple in die Vektoren f und n 
%          erfolgt durch: f=(cos(rak)*cos(str)+cos(dip)*sin(rak)*sin(str),
%                            cos(rak)*sin(str)-cos(dip)*sin(rak)*cos(str),
%                            -sin(rak)*sin(dip)); 
%                         n=(sin(dip)*sin(str),-sin(dip)*cos(str),cos(dip))
%          [siehe hierzu auch Aki/Richards S.106 ff]                   

f=[cos(rak)*cos(str)+cos(dip)*sin(rak)*sin(str),...
      cos(rak)*sin(str)-cos(dip)*sin(rak)*cos(str),...
      -sin(rak)*sin(dip)] 
n=[sin(dip)*sin(str),-sin(dip)*cos(str),cos(dip)]
mom=n'*f+f'*n
ball=zeros(120);
for ix=1:120
   for iy=1:120
      xx=(ix-60)/60;
      yy=(iy-60)/60;
      rq=xx^2+yy^2;
      if rq<1 
         zz=sqrt(1-rq);
         r=[xx;yy;zz;];
         ball(ix,iy)=r'*mom*r;
      end
   end
end
figure(1)
load bal
colormap(bal)
scal=32/max(max(abs(ball)));
ball=max(1,round(ball*scal+32));
image(ball)
axis equal