%close all
%clear all
%----------------------------------------------------------------------------------------------
  n=60;    
  CMT_file = 'jgrb52507-sup-0002-data_s2.txt';         % original USGS FFM data     

  [No Tinit Trise D0 lonR latR dR L0 W0 strike1 dip1 rake1]= textread(CMT_file,...    
   '%d %f %f %f %f %f %f %f %f %f %f %f',n,'headerlines',1);  

  lat2 = NaN*ones(size(No));
  lon2 = NaN*ones(size(No));
  
  zmax   = 100;
  z_plot = zmax*ones(size(D0)); 

FIG1 = figure('Position',[150 50 500 550])

  for i=1:n   
    x1 = lonR(i);    %(ºW) lonitude of Reference Point 
    y1 = latR(i);    %(ºN) latitude of Reference Point 
    L = L0(i);    %(ºW) lonitude of Reference Point 
    W = W0(i);    %(ºN) latitude of Reference Point     
    strike = strike1(i); %(º)  strike angle
    dip = dip1(i);     %(º)  dip angle
  %------------------------------------------------------------------------------
  %
  %   Reference Point
  %         c1_____________c4____\ strike
  %          /            /      / 
  %       W /            /
  %        /____________/
  %      c2      L      c3 
  %
  %       FIG. 1 Fault Area

    %c1 Reference Point
    c1(1) = x1;
    c1(2) = y1;
    %c1 = [x1,y1];

    %c2 : lower left in FIG. 1
    [c2(2),c2(1)] = reckon(y1,x1,km2deg(W)*cos(dip*pi/180),strike+90);

    %c3 : lower right in FIG. 1
    [c3(2),c3(1)] = reckon(c2(2),c2(1),km2deg(L),strike);

    %c4 : upper right in FIG. 1
    [c4(2),c4(1)] = reckon(y1,x1,km2deg(L),strike);

    [cp(2),cp(1)] = reckon(y1,x1,km2deg(L/2),strike);

    x = [c1(1),c2(1),c3(1),c4(1),c1(1)];
    y = [c1(2),c2(2),c3(2),c4(2),c1(2)];
    z = [ D0(i), D0(i), D0(i),  D0(i), D0(i)];
    C=z;

    %3D PLOT
      hFF1 = fill3(x,y,z,z);
      set(gca,'CLim',[0 40])
      hold on

    % colormap default
      colormap jet

      hF1 = plot3(x,y,z,'k-','linewidth',1);
      hold on

      hF2 = plot3(c1(1),c1(2),z_plot,'o');
      set(hF2,'MarkerSize',5,'MarkerFaceColor','w','MarkerEdgeColor','k') 
      hold on 
  
  end
 
  colorbar
 
  axis([140 145 35 41 0 10000]) 
 
 % SET 3D PLOT orientation
   az_rotate = 0.0;
   el_tilt   = 90.0;   
   view(az_rotate, el_tilt);   
  hold on
  xlabel('longitude'), ylabel('latitude')
  
  ind_Tmax = find(Tinit==max(Tinit))   % find the last subfault to move
  Trup     = max(Tinit)+Trise(ind_Tmax) % total rupture duration (sec)
  