% Material visualization with velocity
clear all
close all

dir1='../../';  % input data path
dir2='../../';  % output figure path

% variables on nodes
name = 'RigSubc';         % comp
name_t = 'RigSubt';       % T

% number
nn1 = 500;%[0 10:20:1200];

% extension of data format
ext = '.txt';

% defining size of the model
xsize = 5050;
ysize = 1100;

% defining zoomed area
xbeg  = 1100;
xend  = 3100;
ybeg  = 0;
yend  = 400;

% Defining zoom: YES > 0; NO < 0
zoom = 1;

% Defining melting area
s1 = 0*ones(1,length(nn1));
s2 = 0*ones(1,length(nn1));
s3 = 0*ones(1,length(nn1));
s4 = 0*ones(1,length(nn1)); 

% main loop - insert numbers of images
for idx=1:length(nn1)
    
    cd (dir1)
    
    fname = [name,num2str(nn1(idx)) ext];
    %read input data
    fid = fopen(fname,'r');
    if(fid<=0)
        fname
        disp('No such file...')
        break;
    end
    A = textscan(fid, '%n');
    fclose(fid);
    Data = A{1};
    Time = Data(1);
    nx = Data(2);
    ny = Data(3);
    Data_vec = Data(4:end);
    
    % define coordinates
    dx = xsize/(nx-1);
    dy = ysize/(ny-1);
    gx = 0:dx:xsize;
    gy = 0:dy:ysize;
    s0 = dx*dy;    % the area of grid (km^2) 
   
    if(zoom>0)
        % find beginning/end index
        % x direction
        indx1 = find(gx==xbeg);
        if(length(indx1)==0)
            indx1 = find(gx>xbeg,1);
        end
        indx2 = find(gx==xend);
        if(length(indx2)==0)
            indx2 = find(gx>xend,1);
        end
        % y direction
        indy1 = find(gy==ybeg);
        if(length(indy1)==0)
            indy1 = find(gy>ybeg,1);
        end
        indy2 = find(gy==yend);
        if(length(indy2)==0)
            indy2 = find(gy>yend,1);
        end
        % zoomed area coordinates
        gx = gx(indx1:indx2);
        gy = gy(indy1:indy2);
    end
    
    % initialize composition matrix
    ColorGrid = ones(ny,nx);
    
% colormap for composition
      C_map=[1 1 1;...                       
        0.50588 0.99608 0.78824;...         
        0.4754 0.251 0.2;...                 
        0.68235 0.34118 0;...             
        0.75294 0.75294 0.75294;...        
        211/255 211/255 211/255;...               
        139/250 140/250 140/250;...                     
        13/250 137/250 97/250;...                     
        13/250 137/250 97/250;...               
        160/255 197/255 215/255;...            
        221/255 232/255 236/255;...         
        120/250 178/250 177/250;...        
        12/250 117/250 134/250;...           
        40/255 121/255 168/255;...           
        0.7 0.12 0.3;...                   
        0.4 0.4 0.0;...                     
        255/255 0/255 0/255;...          
        13/250 137/250 97/250;...        
        255/255 0/255 0/255;...             
        0 0 0;...                           
        0 0 0;...                           
        0 0 0;...                           
        0 0 0;...                           
        1 1 0.31765;...                     
        1 0.90196 0.18824;...               
        0.46667 0.46667 0.23529;...         
        0.50196 0.50196 0;...              
        0.8627 0.4902 0.4902;...         
        0.92549 0.43922 0.99608;...         
        13/250 137/250 97/250;...         
        1 1 0;...                         
        1 1 0;...                         
        0.99216 0.38824 0.30196;...         
        0 0 0;...                           
        0 0 0;...                           
        1 0 0;...         
        1 0 0;...        
        1 0 0;...        
        0.8627 0.4902 0.4902;...        
        255/255 0/255 0/255];           
    
       
    %Fill the colorgrid with the data of opened file
   
    num       = 1;
    ind       = 1;
    while num<length(Data_vec)
        value = Data_vec(num);
        if value==-2
            % Compressed: the next ?? indices are given the color material
            num_colors  =   Data_vec(num+1);
            material    =   Data_vec(num+2);
            ind_vec     =   ind:ind+num_colors-1;
            ind         =   ind_vec(end)+1;
            num         =   num+3;
        elseif value==-1
            ind_vec     =   ind;
            material    =   NaN;
            num=num+1;
            ind=ind+1;
        else
            ind_vec     =   ind;
            material    =   value;
            num=num+1;
            ind=ind+1;
        end
        ColorGrid(ind_vec)  =   material;
    end
    
    % zoom in composition
    if(zoom>0)
        ColorGrid = ColorGrid(indy1:indy2,indx1:indx2);
    end
    
    %save memory
    clear Data_vec
    
    time_all(idx)=Time/1e6;
 

    % -------temperature -------
    if(10)
        format long g
        fname_t     =  [name_t,num2str(nn1(idx)) ext];
        fid     =   fopen (fname_t,'r');
        B       =   textscan (fid, '%n')';
        fclose(fid);
        Data    =   B{1};
        nx_t =   Data(2);
        ny_t =   Data(3);
        Data_vec_t=   Data(4:end);
        Data_vec_t= Data_vec_t-273; % form Kelvin to Celcius
        ColorGrid_t = ones(ny_t,nx_t);
        ColorGrid_t = reshape(Data_vec_t,ny_t,nx_t);
        
        % define coordinates for T
        dx_t = xsize/(nx_t-1);
        dy_t = ysize/(ny_t-1);
        gx_t = 0:dx_t:xsize;
        gy_t = 0:dy_t:ysize;
        
        if(zoom>0)
            % find beginning/end index
            % x direction
            indx1_t = find(gx_t==xbeg);
            if(length(indx1_t)==0)
                indx1_t = find(gx_t>xbeg,1);
            end
            indx2_t = find(gx_t==xend);
            if(length(indx2_t)==0)
                indx2_t = find(gx_t>xend,1);
            end
            % y direction
            indy1_t = find(gy_t==ybeg);
            if(length(indy1_t)==0)
                indy1_t = find(gy_t>ybeg,1);
            end
            indy2_t = find(gy_t==yend);
            if(length(indy2_t)==0)
                indy2_t = find(gy_t>yend,1);
            end
            % zoomed area coordinates
            gx_t = gx_t(indx1_t:indx2_t);
            gy_t = gy_t(indy1_t:indy2_t);
        end
        
        % zoom in T
        if(zoom>0)
            ColorGrid_t = ColorGrid_t(indy1_t:indy2_t,indx1_t:indx2_t);
        end
    end 
           
 
 name2 = 'RigSub';

cmyr=100*365*24*3600;

% read data
icount = 0;
ti = 1;
for idx = 1:length(nn1)
    icount = icount+1
    close all
    
    cd(dir1)
    
    fileprn    =  [name2,num2str(nn1(idx)),'.prn'];
    % Opening data file
    fdata=fopen(fileprn,'rb');
    % Read sizes of variables
    A=fread(fdata,4,'uchar');
    % Read model parameters
    % Grid resolution
    xnumx=fread(fdata,1,'int64');
    ynumy=fread(fdata,1,'int64');
    % Markers per cell
    mnumx=fread(fdata,1,'int64');
    mnumy=fread(fdata,1,'int64');
    % Number of markers
    marknum=fread(fdata,1,'int64');
    % Model sizes
    xsize=fread(fdata,1,'float64');
    ysize=fread(fdata,1,'float64');
    % Pressure value
    pinit=fread(fdata,5,'float64');
    % Gravity
    GXKOEF=fread(fdata,1,'float64');
    GYKOEF=fread(fdata,1,'float64');
    % Number of rocks
    rocknum=fread(fdata,1,'int32');
    % Number of Boundary conditions
    bondnum=fread(fdata,1,'int64');
    % Stage,time
    n1=fread(fdata,1,'int32');
    timesum=fread(fdata,1,'float64');
    if(timesum>1e3)
        time=sprintf('%3.3f',timesum/1e6);
        time=str2num(time);
    else
        time=timesum/1e6;
    end
    time_all(icount) = time; % Ma
    
    % Skip general information
    curpos0 = 4*1;
    % skip grid parameters
    curpos1 = 12*8 + 2*8 + 4 + 8 + 4 + 8;
    % Skip rock informations
    curpos2 = rocknum*(1*4 + 24*8);
    % Skip nodes information
    curpos3 = xnumx*ynumy*(3*4 + 3*8 + 16*4 + 1*8 + 5*4);        %change
    % skip gridline positions
    curpos4 = (xnumx + ynumy)*4;
    % skip boundary conditions equations
    curpos5 = (bondnum-1)*(4*4 + 3*8);
    % skip marker infos
    curpos6 = (marknum+1)*(9*4+1);
    
    % total value should equals to the real prn size (in byte)
    curpos = curpos0 + curpos1 + curpos2 + curpos3 + curpos4 + curpos5 + curpos6;
    
    % initialitialise parameter
    sum_nu = 0; % mean effective viscosity
    sum_eii = 0;
    
    
    %------- ------------load gridlines ---------------------------
    curpos_grid = curpos0 + curpos1 + curpos2 + curpos3;
    fseek(fdata,curpos_grid,'bof')
    gx_gy = fread(fdata,xnumx+ynumy,'float32');
    gx2 = gx_gy(1:xnumx)/1000;
    gy2 = gx_gy(xnumx+1:end)/1000;
    
    % ------------ fast: write data into a 2d matrix -------------
    % find the starting position
    % sxx
    curpos_sxx = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 3*4;
    fseek(fdata,curpos_sxx,'bof');
    % only read sxx, skip useless information
    skip1 = 3*4 + 3*8 + 16*4 + 1*8 + 4*4;                         %change
    % 2d matrix
    sxx = fread(fdata, [ynumy,xnumx],'float32',skip1);
    
    % syy
    curpos_syy = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 4*4;
    fseek(fdata,curpos_syy,'bof');
    % only read sxx, skip useless information
    skip_syy = skip1;
    % 2d matrix
    syy = fread(fdata, [ynumy,xnumx],'float32',skip_syy);
    
    % sxy
    curpos_sxy = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 5*4;
    fseek(fdata,curpos_sxy,'bof');
    skip2 = skip1;
    sxy = fread(fdata, [ynumy,xnumx],'float32',skip2);
    
    % exx
    curpos_exx = curpos0 + curpos1 + curpos2 + 3*4+3*8;
    fseek(fdata,curpos_exx,'bof');
    skip3 = skip1;
    exx = fread(fdata, [ynumy,xnumx],'float32',skip3);
    
    % eyy
    curpos_eyy = curpos0 + curpos1 + curpos2 + 3*4+3*8 +4;
    fseek(fdata,curpos_eyy,'bof');
    skip3 = skip1;
    eyy = fread(fdata, [ynumy,xnumx],'float32',skip3);
    
    % exy
    curpos_exy = curpos0 + curpos1 + curpos2 +  3*4+3*8 + 2*4;
    fseek(fdata,curpos_exy,'bof');
    skip4 = skip1;
    exy = fread(fdata, [ynumy,xnumx],'float32',skip4);
    
    % vx
    curpos_vx = curpos0 + curpos1 + curpos2 + 4;
    fseek(fdata,curpos_vx,'bof');
    skip3 = skip1;
    vx = fread(fdata, [ynumy,xnumx],'float32',skip3);
    vx = vx*cmyr;
    
    % vy
    curpos_vy = curpos0 + curpos1 + curpos2 + 2*4;
    fseek(fdata,curpos_vy,'bof');
    skip3 = skip1;
    vy = fread(fdata, [ynumy,xnumx],'float32',skip3);
    vy = vy*cmyr;
    
    % RMS velocity
    rms_vel(icount) = ((sum(sum(vx.^2)) + sum(sum(vy.^2)))/(2*xnumx*ynumy))^0.5;
    rms_vx(icount)  = ((sum(sum(vx.^2)))/(xnumx*ynumy))^0.5;
    rms_vy(icount)  = ((sum(sum(vy.^2)))/(xnumx*ynumy))^0.5;
    
    % shear heating
    Hs = 2*sxx.*exx + 2*sxy.*exy;
    aa = find(Hs<1e-30);
    Hs(aa) = 1e-30;
    
    % density
    curpos_dens = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 6*4;
    fseek(fdata,curpos_dens,'bof');
    dens = fread(fdata, [ynumy,xnumx],'float32',skip1);
    
    % second invariant
    sii = (sxx.^2 + sxy.^2).^(1/2);
    eii = (exx.^2*0.5 + eyy.^2*0.5 + exy.^2).^(1/2);
    eii(find(eii<=0)) = 1e-30;
    
    % viscosity
    curpos_nu = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 7*4;
    fseek(fdata,curpos_nu,'bof');
    skip5 = skip1;
    nu = fread(fdata, [ynumy,xnumx],'float32',skip5);
    
    % pr(pressure)
    curpos_pr = curpos0 + curpos1 + curpos2;
    fseek(fdata,curpos_pr,'bof');
    %only read sxx, skip useless information
    skip3 = skip1;%跳转到下一个节点该数据的位置进行读取
    %2d matrix
    pr = fread(fdata,[ynumy,xnumx],'float32',skip3);
    
    % temperature
    curpos_T = curpos0 + curpos1 + curpos2 + 3*4+3*8+ 15*4;
    fseek(fdata,curpos_T,'bof');
    T = fread(fdata, [ynumy,xnumx],'float32',skip1);
    T=T-273;
    % mext
    curpos_mext = curpos_T + 8+4*4;        %change
    fseek(fdata,curpos_mext,'bof');                         %change
    mext = fread(fdata, [ynumy,xnumx],'float32',skip1);    %change
     
    % meltfrac
    curpos_meltfrac = curpos_T +8+5*4;        %change
    fseek(fdata,curpos_meltfrac,'bof');                         %change
    meltfract = fread(fdata, [ynumy,xnumx],'float32',skip1);    %change
    
    
    % mean viscosity
    % entire box
    mean_nu_box(idx) = mean(mean(nu));
    % mantle lithosphere
    nodecount1 = find(T<=1300);
    nodecount2 = find(dens>3000);
    nodecount3 = intersect(nodecount1,nodecount2);
    mean_nu_lith(idx) = mean(nu(nodecount3));
    % asthenosphere
    nodecount4 = find(T>1300);
    mean_nu_asth(idx) = mean(nu(nodecount4));
    

    % ----------------------- Plot -------------------------
    lw = 1.2;
    ft = 12;

    
    % ------------------------ comp plot ------------------------
    if(10)
        cd (dir2)
        figure(1)
        step=1;
        colormap(C_map);
        pcolor(gx,gy,ColorGrid);
        shading interp
       if(10)
        hold on
           [C,h] = contour(gx_t(1:step:end),gy_t(1:step:end),ColorGrid_t(1:step:end,1:step:end),...
               [300:300:1200],'color','[0.5 0.5 0.5]','LineWidth',0.2);
            text_handle = clabel(C,h,'color','k','labelspacing',10000,'fontsize',6);
        hold off
       end
       
       if(10)
           step1=15;
            hold on;
             yyid=190;
            scale= 10;
                 h1=quiver(gx2(70:step1:430), gy2(20:step1:yyid), vx(20:step1:yyid,70:step1:430), vy(20:step1:yyid,70:step1:430),0,'k','linewidth',0.5);
                hu=get(h1,'UData'); hv=get(h1,'VData');
                 set(h1,'UData',scale*hu,'VData',scale*hv);
             hold off;       
       end
         caxis ([0 38])
         axis ij equal tight
              xlim ([1100 3100])
              ylim ([0 400])
              set(gca,'xtick',1200:200:3000);
          title(['Time = ', num2str(Time/1e6),' Ma'],'fontsize',ft)
          set(gca,'position',[0 0 1 1])
        if zoom >= 0
            filename    =  ['Comp_z_prn',num2str(nn1(idx)),'.png'];
             set(gca,'layer','top','position',[0.07 0.08 0.9 0.9],'linewidth',0.01)  
        else
            filename    =  ['Comp_prn',num2str(nn1(idx))];
            set(gca,'layer','top','linewidth',lw,'fontsize',ft)
        end
         print ('-dpng', '-r800',filename);
    end
      
end
end



