output_plot_blocks1

PURPOSE ^

Compile, plot and save the results after optimization

SYNOPSIS ^

function output_plot_blocks1(Blocks120, Result, numOfPeriods, numBlocks,normC, xdim, ydim, zdim, v, p, q,miningGoal, processingGoal, dykemGoal,miningd1, processingd2, dykemd3, dykemd4)

DESCRIPTION ^

 Compile, plot and save the results after optimization 
--------------------------------------------------------------------------
 Date: July 30, 2010, ver01
 By: Eugene Ben-Awuah
--------------------------------------------------------------------------

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % Compile, plot and save the results after optimization
0002 %--------------------------------------------------------------------------
0003 % Date: July 30, 2010, ver01
0004 % By: Eugene Ben-Awuah
0005 %--------------------------------------------------------------------------
0006 
0007 function output_plot_blocks1(Blocks120, Result, numOfPeriods, numBlocks,...
0008                           normC, xdim, ydim, zdim, v, p, q,...
0009                           miningGoal, processingGoal, dykemGoal,...
0010                           miningd1, processingd2, dykemd3, dykemd4) 
0011                       
0012 % steps that will take to plot
0013 step = 3;
0014  
0015 % xdim = dimension of block in the X direction, Easting;
0016 % ydim = dimension of block in the Y direction, Northing;
0017 % zdim = dimension of block in the Z direction, Elevation;
0018  
0019 % algorithm used by sovler
0020 solverAlgorithm = Result.SolverAlgorithm 
0021  
0022 % CPU time used in seconds
0023 cpuTime = Result.CPUtime 
0024  
0025 % number of major iterations
0026 % numOfIterations = Result.iter
0027  
0028 % Number of constraint evaluations needed.
0029 constraintEvalualtions = Result.ConstrEv
0030  
0031 % discounted cash flow at the optimum schedule
0032 display('Discounted Cash Flow in Million Dollars')
0033 DCF = cur2str((Result.f_k*-1* normC)/1000000,2)
0034  
0035 %--------------------------------------------------------------------------
0036 % binSchedule [numBlocks, numOfPeriods]= binary integer matrix holding the
0037 % schedule
0038 % each column contains the corresponding binary schedule of that period.
0039 % e.g. second column holds the second period schedule
0040  
0041 % b(K*T) = Result.x_k( 1: K*T )= holds the binary integer variable
0042 % y(K*T) = Result.x_k( [end-(K*T)+1] : end )= holds the continuous cut
0043 %                                             variable
0044 % x(N*T) = Result.x_k( [(K*T)+1] : [end-(K*T)] )= holds the continuous block
0045 %                                                 variable
0046  
0047 x  = Result.x_k( (numBlocks*numOfPeriods+1) : (2*(numBlocks*numOfPeriods)));
0048 z  = Result.x_k( ((2*numBlocks*numOfPeriods)+1) : (3*(numBlocks*numOfPeriods))); 
0049 y  = Result.x_k( ((3*numBlocks*numOfPeriods)+1) : (4*(numBlocks*numOfPeriods)));
0050 d1 = Result.x_k( ((4*numBlocks*numOfPeriods)+1) : ((4*(numBlocks*numOfPeriods)+numOfPeriods)));
0051 d2 = Result.x_k( ((4*(numBlocks*numOfPeriods)+numOfPeriods)+1) : ((4*(numBlocks*numOfPeriods)+(2*numOfPeriods))));
0052 d3 = Result.x_k( ((4*(numBlocks*numOfPeriods)+(2*numOfPeriods))+1) : ((4*(numBlocks*numOfPeriods)+(3*numOfPeriods))));
0053 d4 = Result.x_k( ((4*(numBlocks*numOfPeriods)+(3*numOfPeriods))+1) : ((4*(numBlocks*numOfPeriods)+(4*numOfPeriods))));
0054 
0055 binOre_Schedule = reshape(x, numBlocks, numOfPeriods);
0056 binDyke_Schedule = reshape(z, numBlocks, numOfPeriods); 
0057 binMining_Schedule = reshape(y, numBlocks, numOfPeriods);
0058 d1 = d1'; d2 = d2'; d3 = d3'; d4 = d4';
0059 binWaste_Schedule = binMining_Schedule - binOre_Schedule - binDyke_Schedule;
0060 % for iPeriods = 1: numOfPeriods
0061 %   binSchedule(:,iPeriods) = Result.x_k((numCuts * (iPeriods))+ 1:...
0062 %               (numCuts * (iPeriods-1))+ numBlocks ,1);
0063 % end
0064  
0065 %[m n]=size(binOre_Schedule);
0066 %binWasteSchedule = zeros(m, n);
0067 %binOreSchedule = zeros(m, n);
0068 %binDykeSchedule = zeros(m, n);
0069  
0070 % constructing matricies with the same format of binSchedule used for plots
0071 for iPeriods = 1: numOfPeriods
0072     for iBlocks = 1: numBlocks
0073     blockTonnage(iBlocks,iPeriods) = Blocks120(iBlocks).BlockTonnage;
0074     EBV(iBlocks,iPeriods) = Blocks120(iBlocks).EBV;
0075  
0076     %if (EBV(iBlocks,iPeriods)<=0 && binSchedule(iBlocks,iPeriods)~= 0)
0077         % binWasteSchedule(iBlocks,iPeriods) = 1;
0078         %binWasteSchedule(iBlocks,iPeriods) = binSchedule(iBlocks,iPeriods);
0079      %elseif (EBV(iBlocks,iPeriods)> 0 && binSchedule(iBlocks,iPeriods)~= 0)
0080          %binOreSchedule(iBlocks,iPeriods) = 1;
0081          %binOreSchedule(iBlocks,iPeriods) = binSchedule(iBlocks,iPeriods);
0082      %end
0083  
0084      blockBitumenOre(iBlocks,iPeriods) = Blocks120(iBlocks).blockBitumenOre;
0085      blockFinesOre(iBlocks,iPeriods) = Blocks120(iBlocks).blockFinesOre;
0086      %P(iBlocks,iPeriods) = Blocks(iBlocks).P;
0087      bitumenGradeOre(iBlocks,iPeriods) = Blocks120(iBlocks).bitumenGradeOre;
0088      finesGradeOre(iBlocks,iPeriods) = Blocks120(iBlocks).finesGradeOre;
0089      
0090      blockFinesDyke(iBlocks,iPeriods) = Blocks120(iBlocks).blockFinesDyke;
0091      %P(iBlocks,iPeriods) = Blocks(iBlocks).P;
0092      finesGradeDyke(iBlocks,iPeriods) = Blocks120(iBlocks).finesGradeDyke;
0093      %gradeP(iBlocks,iPeriods) = Blocks(iBlocks).gradeP;
0094      OreTonnes(iBlocks,iPeriods) = Blocks120(iBlocks).oreTonnage;
0095      DykeTonnes(iBlocks,iPeriods) = Blocks120(iBlocks).dykeTonnage;
0096      WasteTonnes(iBlocks,iPeriods) = Blocks120(iBlocks).wasteTonnage;
0097      OreValue(iBlocks,iPeriods) = Blocks120(iBlocks).oreValueCutoff;
0098      DykeCost(iBlocks,iPeriods) = Blocks120(iBlocks).dykeCostCutoff;
0099      WasteCosts(iBlocks,iPeriods) = Blocks120(iBlocks).miningCost;
0100      
0101   end
0102  end
0103  
0104      tonnageMined = binMining_Schedule .* blockTonnage;
0105      tonnageMined1 = sum(sum(tonnageMined));
0106      
0107      % NPV at the optimum schedule
0108      display('NPV in Million Dollars')
0109      NPV = (binOre_Schedule.*v) + (binDyke_Schedule.*p)+ (binMining_Schedule.*q);
0110      NPV = sum(sum(NPV));
0111      NPV = cur2str((NPV)/1000000,2)
0112      TotalDykeCost = sum(sum(binDyke_Schedule.*p))
0113      %cashflow1 = sum(cashflow)
0114      cashFlow = binMining_Schedule.* EBV; % + binDyke_Schedule.* DykeCost;
0115      %Modify the EBV values to reflect the material segregation implemented
0116      %using cutoff analysis (ore,dyke and waste material)
0117      %Overall_cashFlow = binOre_Schedule.* EBV + binDyke_Schedule.* DykeCost...
0118      %                   + binWaste_Schedule.* EBV;
0119      
0120      % ore send to mill
0121      % productTonnage = binSchedule .* MWT;
0122      oreTonnage = binOre_Schedule .* OreTonnes;
0123      %%oreTonnage = binMining_Schedule .* OreTonnes;
0124      oreTonnage1 = sum(sum(oreTonnage));
0125      finesOreTonnage = binOre_Schedule .* blockFinesOre;
0126      %%finesOreTonnage = binMining_Schedule .* blockFines;
0127      
0128      %%dykeTonnage = binMining_Schedule .* DykeTonnes;
0129      dykeTonnage = binDyke_Schedule .* DykeTonnes;
0130      dykeTonnage1 = sum(sum(dykeTonnage));
0131      %%dykeTonnage2 = binDyke_Schedule .* blockTonnage
0132      finesDykeTonnage = binDyke_Schedule .* blockFinesDyke;
0133      %%finesDykeTonnage = binMining_Schedule .* blockFines;
0134      
0135      %%wasteTonnage = binWaste_Schedule .* WasteTonnes;
0136      %%wasteTonnage2 = sum(sum(wasteTonnage2))
0137      %%%wasteTonnage = binMining_Schedule .* WasteTonnes;
0138      %wasteTonnage = binWaste_Schedule .* blockTonnage;
0139      %%wasteTonnage = sum(sum(wasteTonnage))
0140      %!!!Need to investigate why there are differences between wasteTonnage
0141      %and wasteTonnage2
0142      %%This reconciles with wasteTonnage
0143      wasteTonnage = tonnageMined - oreTonnage - dykeTonnage; 
0144      wasteTonnage1 = sum(sum(wasteTonnage));
0145      display('Scheduled Materials  in Million Tonnes')
0146      tonnageMined1 = tonnageMined1/1000000
0147      oreTonnage1 = oreTonnage1/1000000
0148      dykeTonnage1 = dykeTonnage1/1000000
0149      wasteTonnage1 = wasteTonnage1/1000000
0150      %sTonnage = binSchedule .* S;
0151      %pTonnage = binSchedule .* P;
0152      
0153  %   oreGrade = binSchedule .*gradeMWT;
0154  %   sGrade = binSchedule .* gradeS;
0155  %   pGrade = binSchedule .* gradeP;
0156   
0157      %oreGrade = gradeMWT;
0158      %sGrade = gradeS;
0159      %pGrade = gradeP;
0160  
0161  % Plotting grades, etc-----------------------------------------------------
0162  % createBar(data, xLabel, yLabel, Title)
0163  
0164  createBar(sum(tonnageMined)/1000000,...
0165                   'Period',...
0166                   'Tonnage Mined (MTonne)',...
0167                   'Total Tonnage Mined');
0168 hold on; 
0169 mG=plot(miningGoal);
0170 set(mG,'Color','red','LineWidth',2)
0171 hold off;    
0172 
0173  createBar(sum(oreTonnage)/1000000,...
0174                  'Period',...
0175                  'Ore Tonnage (MTonne)',...
0176                  'Ore Tonnage');
0177 hold on; 
0178 pG=plot(processingGoal);
0179 set(pG,'Color','red','LineWidth',2);
0180 hold off;    
0181 
0182 createBar(sum(dykeTonnage)/1000000,...
0183                  'Period',...
0184                  'DykeM Tonnage (MTonne)',...
0185                  'Dyke Material Tonnage');
0186 hold on; 
0187 dG=plot(dykemGoal);
0188 set(dG,'Color','red','LineWidth',2);
0189 hold off;  
0190 
0191 createBar(sum(wasteTonnage)/1000000,...
0192                  'Period',...
0193                  'Waste Tonnage (MTonne)',...
0194                  'Waste Tonnage');
0195                                                            
0196  Y = [sum(oreTonnage)/1000000; sum(dykeTonnage)/1000000;... 
0197       sum(wasteTonnage)/1000000]';
0198  createBar(Y,...
0199              'Period',...
0200              'Ore, dyke and waste tonnage (MTonne)',...
0201              'Tonnage of Ore, Dyke and Waste');
0202                
0203  createBar((d1)/1000000,...
0204                  'Period',...
0205                  'Miningd1 Tonnage (MTonne)',...
0206                  'Mining negative deviation tonnage');
0207 hold on; 
0208 md1=plot(miningd1);
0209 set(md1,'Color','red','LineWidth',2);
0210 hold off;  
0211 
0212 %Create errorbars for mining, processing and dyke material using the
0213 %deviational variables used in the optimization
0214 ebperiod = 1:numOfPeriods;
0215 ebU = zeros(1, numOfPeriods);
0216 ebL = d1/1000000;
0217 errorbar(ebperiod, sum(tonnageMined)/1000000, ebU, ebL)
0218 xlabel('Period');
0219 ylabel('Tonnage Mined (MTonne)');
0220 title('Tonnage mined and deviation from target'); 
0221 
0222 createBar((d2)/1000000,...
0223                  'Period',...
0224                  'Processingd2 Tonnage (MTonne)',...
0225                  'Processing negative deviation tonnage');
0226 hold on; 
0227 pd2=plot(processingd2);
0228 set(pd2,'Color','red','LineWidth',2);
0229 hold off;  
0230 
0231 ebL = d2/1000000;
0232 errorbar(ebperiod, sum(oreTonnage)/1000000, ebU, ebL);
0233 xlabel('Period');
0234 ylabel('Ore Tonnage (MTonne)');
0235 title('Ore Tonnage and deviation from target');
0236 
0237 createBar((d3)/1000000,...
0238                  'Period',...
0239                  'Dykemd3 Tonnage (MTonne)',...
0240                  'Dyke material negative deviation tonnage');
0241 hold on; 
0242 dm3=plot(dykemd3);
0243 set(dm3,'Color','red','LineWidth',2);
0244 hold off;  
0245 
0246 createBar((d4)/1000000,...
0247                  'Period',...
0248                  'Dykemd4 Tonnage (MTonne)',...
0249                  'Dyke material positive deviation tonnage');
0250 hold on; 
0251 dm4=plot(dykemd4);
0252 set(dm4,'Color','red','LineWidth',2);
0253 hold off;  
0254 
0255 ebL = d3/1000000;
0256 ebU = d4/1000000;
0257 errorbar(ebperiod, sum(dykeTonnage)/1000000, ebU, ebL);
0258 xlabel('Period');
0259 ylabel('Dyke material tonnage (MTonne)');
0260 title('Dyke material tonnage and deviation from target');
0261 
0262 createBar(sum(cashFlow)/1000000,...
0263                   'Period',...
0264                   'cash flow ($M)',...
0265                   'Cash Flow');
0266               
0267   %createBar(sum(Overall_cashFlow)/1000000,...
0268   %                'Period',...
0269   %                'overall cash flow ($M)',...
0270   %                'Overall Cash Flow');
0271  % %----------------------------------------------------------------------
0272  % UNCOMMENT this section to plot high resoultion Figures
0273  % Plotting the same figures as above with high resolution as TIF files
0274  %
0275  % createBar(sum(tonnageMined)/1000000,...
0276  %                  'Period',...
0277  %                  'Tonnage Mined(MTonne)',...
0278  %                  'Total Tonnage Mined');
0279  %
0280  % print -dtiff -r150 im.tiff
0281  % im = imread( 'im.tiff', 'tiff' );
0282  % [im_hatch,colorlist] = applyhatch_pluscolor(im,'\-x.',1,[],[],100,3,3);
0283  % imwrite(im_hatch,'im_hatch.tiff','tiff')
0284  %
0285  % createBar(sum(oreTonnage)/1000000,...
0286  %                 'Period',...
0287  %                 'Ore Tonnage (MTonne)',...
0288  %                 'Ore Tonnage');
0289  %
0290  % print -dtiff -r150 im.tiff
0291  % im = imread( 'im.tiff', 'tiff' );
0292  % [im_hatch,colorlist] = applyhatch_pluscolor(im,'\-x.',1,[],[],100,3,3);
0293  % imwrite(im_hatch,'im_hatch.tiff','tiff')
0294  %
0295  %  createBar(sum(wasteTonnage)/1000000,...
0296  %                 'Period',...
0297  %                 'Waste Tonnage (MTonne)',...
0298  %                 'Waste Tonnage');
0299  %
0300  % print -dtiff -r150 im.tiff
0301  % im = imread( 'im.tiff', 'tiff' );
0302  % [im_hatch,colorlist] = applyhatch_pluscolor(im,'\-x.',1,[],[],100,3,3);
0303  % imwrite(im_hatch,'im_hatch.tiff','tiff')
0304  %
0305  % Y = [sum(oreTonnage)/1000000; sum(wasteTonnage)/1000000]';
0306  % createBar(Y,...
0307  %             'Period',...
0308  %             'Ore and waste tonnage (MTonne)',...
0309  %             'Tonnage of Ore and Waste');
0310  %
0311  % print -dtiff -r150 im.tiff
0312  % im = imread( 'im.tiff', 'tiff' );
0313  % [im_hatch,colorlist] = applyhatch_pluscolor(im,'\-x.',1,[],[],100,3,3);
0314  % imwrite(im_hatch,'im_hatch.tiff','tiff')
0315  %
0316  % createBar(sum(cashFlow)/1000000,...
0317  %                  'Period',...
0318  %                  'cash flow (M$)',...
0319  %                  'Cash Flow');
0320  %
0321  % print -dtiff -r150 im.tiff
0322  % im = imread( 'im.tiff', 'tiff' );
0323  % [im_hatch,colorlist] = applyhatch_pluscolor(im,'\-x.',1,[],[],100,3,3);
0324  % imwrite(im_hatch,'im_hatch.tiff','tiff')
0325  % %------------------------------------------------------------------------
0326   
0327   for iPeriods = 1: numOfPeriods
0328         
0329         % SumProduct of ore grade and tonnage to calculate average grade
0330         sumProductGradeOreTonnage(iPeriods) = ...
0331                          sum(bitumenGradeOre(:,iPeriods).* oreTonnage(:,iPeriods));
0332         % sum of ore tonnage
0333         sumOreTonnage(iPeriods) = sum(oreTonnage(:,iPeriods));
0334         
0335         sumProductFinesOreTonnage(iPeriods) = ...
0336                          sum(finesGradeOre(:,iPeriods).* finesOreTonnage(:,iPeriods));
0337         sumFinesOreTonnage(iPeriods) = sum(finesOreTonnage(:,iPeriods));
0338         
0339         sumProductFinesDykeTonnage(iPeriods) = ...
0340                          sum(finesGradeDyke(:,iPeriods).* finesDykeTonnage(:,iPeriods));
0341         sumFinesDykeTonnage(iPeriods) = sum(finesDykeTonnage(:,iPeriods));
0342         
0343     end 
0344  %--------------------------------------------------------------------------
0345  % createPlot(xData, yData, xLabel, yLabel, Title)
0346  
0347   createPlot( 1:numOfPeriods,...
0348               (sumProductGradeOreTonnage./sumOreTonnage)*100,...
0349               'Period',... 
0350               ' bit% ',...
0351               'Average ore grade bit% per period');
0352  
0353   createPlot( 1:numOfPeriods,...
0354               (sumProductFinesOreTonnage./sumFinesOreTonnage)*100,...
0355               'Period',...
0356               'fines%',...
0357               'Average ore grade fines% per period');
0358               
0359   createPlot( 1:numOfPeriods,...
0360               (sumProductFinesDykeTonnage./sumFinesDykeTonnage)*100,...
0361               'Period',...
0362               'fines% ',...
0363               'Average dyke grade fines per period');             
0364               
0365  % struct Graphs holds the attributes per year------------------------------
0366      Graphs = struct(...
0367          'oreTonnagePerYr',sum(oreTonnage)/1000000,...
0368          'dykeTonnagePerYr',sum(dykeTonnage)/1000000,...
0369          'wasteTonnagePerYr' , sum(wasteTonnage)/1000000,...
0370          'cashFlowPerYr', sum(cashFlow)/1000000,...
0371          'bitOreGradePerYr', (sumProductGradeOreTonnage./sumOreTonnage)*100,...
0372          'finesOreGradePerYr',   (sumProductFinesOreTonnage./sumFinesOreTonnage)*100,...
0373          'finesDykeGradePerYr',  (sumProductFinesDykeTonnage./sumFinesDykeTonnage)*100,...
0374          'tonnageMinedPerYr',(sum(tonnageMined)/1000000),...
0375          'miningNegativeDeviationPerYr',(d1/1000000),...
0376          'processingNegativeDeviationPerYr',(d2/1000000),...
0377          'dykemNegativeDeviationPerYr',(d3/1000000),...
0378          'dykemPostiveDeviationPerYr',(d4/1000000));
0379          
0380      save graphs Graphs;
0381     
0382  for iBlocks = 1: numBlocks
0383         
0384   % adds the array Blocks.periods(1:2,:) to the Blocks struct
0385   % Blocks.periods(1 ,:) = periods that each block is extracted in
0386   % Blocks.periods(2 ,:) = portion of the block that is extracted in the
0387   % corrosponding period.
0388   
0389   % returns the periods
0390   Blocks120(iBlocks).periods(1, :) = find (binMining_Schedule(iBlocks, :));
0391   
0392    % returns the portion of the block extracted in each period
0393   Blocks120(iBlocks).periods(2, :) = binMining_Schedule(iBlocks, ...
0394                                              Blocks120(iBlocks).periods(1, :));
0395  end % end for iBlocks
0396      
0397  % save the Blocks with their respective extraction periods.
0398  save scheduledBlocks Blocks120      
0399  
0400  %--------------------------------------------------------------------------
0401  % finding the maximum and minimum block indices to plot the schedule
0402  x = zeros(numBlocks, 1); y = zeros(numBlocks, 1); z = zeros(numBlocks, 1);
0403  
0404  for iBlocks = 1: numBlocks
0405      x(iBlocks) = Blocks120(iBlocks).XI;
0406      y(iBlocks) = Blocks120(iBlocks).YI;
0407      z(iBlocks) = Blocks120(iBlocks).ZI;
0408  end % end for iBlocks
0409  
0410  % max_XI holds the maximum index of blocks in X direction
0411  % max_XI_index holds the maximimum general index (iBlocks = 1:numBlocks)
0412  % of block max_XI
0413  [max_XI, max_XI_index] = max(x); [min_XI, min_XI_index] = min(x);
0414  [max_YI, max_YI_index] = max(y); [min_YI ,min_YI_index] = min(y);
0415  [max_ZI, max_ZI_index] = max(z); [min_ZI, min_ZI_index] = min(z);
0416  
0417  % coordinates of the origin of the blocks that will be plotted
0418  xCoorMin = Blocks120(min_XI_index).X;
0419  yCoorMin = Blocks120(min_YI_index).Y;
0420  zCoorMin = Blocks120(min_ZI_index).Z;
0421  
0422  xCoorMax = Blocks120(max_XI_index).X;
0423  yCoorMax = Blocks120(max_YI_index).Y;
0424  zCoorMax = Blocks120(max_ZI_index).Z;
0425  %--------------------------------------------------------------------------
0426  % switching the Matlab matrix coordinates to match the whittle coordinates
0427  % check documentation for the reason
0428   max_matrix_XI = max_YI;
0429   max_matrix_YI = max_XI;
0430   max_matrix_ZI = max_ZI;
0431   
0432   min_matrix_XI = min_YI;
0433   min_matrix_YI = min_XI;
0434   min_matrix_ZI = min_ZI;
0435     
0436   blockModel = (numOfPeriods + 5) * ones( max_matrix_XI,...
0437                                           max_matrix_YI,...     
0438                                           max_matrix_ZI );     
0439  for iBlocks = 1: numBlocks
0440      
0441      % see documentation for this section
0442      % XI location based on Whittle coordinate
0443      % matrix_XI coordinates
0444      
0445      % these are for checking delete later on
0446      % Blocks(iBlocks).XI;
0447      % Blocks(iBlocks).YI;
0448      % Blocks(iBlocks).ZI;
0449      
0450        matrix_XI = max_matrix_XI - Blocks120(iBlocks).YI + 1;
0451        matrix_YI = Blocks120(iBlocks).XI;
0452        matrix_ZI = Blocks120(iBlocks).ZI;
0453      
0454      % this is for checking delete later on
0455      %  period = Blocks(iBlocks).periods(1,1)
0456     
0457      % For the purpose of visualizing the blocks a decision must be made
0458      % for the blocks that are extracted in different periods. here we are
0459      % consdidering the period that has the largest fraction for
0460      % visualizing this can cause some problems in the order of the blocks
0461      % that are shown! be aware and make sure it is working fine.
0462      
0463      %This checks and assigns period and fraction mined to blocks that are not
0464      %mined at all. These blocks are assigned empty matrices by the code
0465      %and creates errors during ploting
0466      if isempty(Blocks120(iBlocks).periods)==1;
0467              Blocks120(iBlocks).periods=[(numOfPeriods + 5);1];
0468      end
0469      
0470       [maxVal index] = max(Blocks120(iBlocks).periods(2,:));
0471      
0472       blockModel(matrix_XI, matrix_YI, matrix_ZI) = ...
0473                                            Blocks120(iBlocks).periods(1,index);
0474      
0475  %     % OR the first period that a portion of a block is extracted
0476  %     blockModel(matrix_XI, matrix_YI, matrix_ZI) = ...
0477  %                                            Blocks(iBlocks).periods(1,end);
0478  
0479      % save the extraction period in Period to export to the excel file
0480      Period(iBlocks) = Blocks120(iBlocks).periods(1,1);
0481         
0482  end % end for iBlocks
0483  
0484  % transpose
0485  % Period(1:numBlocks, 1)= holds the period that block i has been last
0486  % extracted
0487  Period = Period';
0488  save Period Period;
0489  
0490  blockModel = blockModel( (max_YI - max_matrix_XI + 1):...
0491                           (max_matrix_XI - min_YI + 1),...
0492                            min_matrix_YI: max_matrix_YI,...
0493                            min_matrix_ZI: max_matrix_ZI );
0494  
0495   % m in the Y direction - Northing
0496   % n in the X direction - Easting
0497   % r in the Z direction - Eleveation
0498   [m n r] = size(blockModel);                       
0499   
0500   % this is how the block model is setup
0501   % blockModel(Y, X, Z)
0502  
0503  % Elevation - Planview
0504  for k = 1:r
0505      temp = blockModel(:,:,k);
0506      temp = squeeze(temp);
0507      figure
0508      % axis([xCoorMin xCoorMax yCoorMin yCoorMax]);
0509      viewmatrix(temp)
0510      xlabel( {'Easting (X-index)'; ['each cell ' num2str(xdim) 'm']},...
0511            'FontName','Times New Roman',...
0512            'FontSize',14);
0513      ylabel( {'Northing (Y-index)'; ['each cell ' num2str(ydim) 'm']},... 
0514            'FontName','Times New Roman',...
0515            'FontSize',14);
0516      title(['Plan View - ' num2str(zCoorMin + (k-1)*zdim) 'm'],...
0517             'FontName','Times New Roman',...
0518             'FontSize',14);
0519      colormap(hot)
0520      set(gca,'DataAspectRatio',[ydim xdim 1])     
0521  end
0522  
0523  % cross sections looking East on the X axis
0524  for j = 1:step:n
0525      temp = blockModel(:,j,:);
0526      temp = squeeze(temp);
0527      temp = rot90(temp);
0528      figure
0529      viewmatrix(temp)
0530      xlabel({'Northing (Y-index)'; ['each cell ' num2str(ydim) 'm']},... 
0531             'FontName','Times New Roman',...
0532             'FontSize',14);
0533      ylabel({'Elevation (Z-index)'; ['each cell ' num2str(zdim) 'm']},... 
0534              'FontName','Times New Roman',...
0535              'FontSize',14);
0536      title(['Cross Section Looking East - '...
0537              num2str(xCoorMin + (j-1)*xdim) 'm'],...
0538             'FontName','Times New Roman',...
0539             'FontSize',14);
0540      colormap(hot)
0541      set(gca,'DataAspectRatio',[zdim ydim 1])
0542  end
0543  
0544  % cross sections looking North on the Y axis
0545  for i = 1:step:m
0546      temp = blockModel(i,:,:);
0547      temp = squeeze(temp);
0548      temp = rot90(temp);
0549      figure
0550      viewmatrix(temp)
0551      xlabel({'Easting (X-index)'; ['each cell ' num2str(xdim) 'm']},...
0552              'FontName','Times New Roman',...
0553              'FontSize',14);
0554      ylabel({'Elevation (Z-index)'; ['each cell ' num2str(zdim) 'm']},... 
0555            'FontName','Times New Roman',...
0556            'FontSize',14);
0557      title(['Cross Section Looking North - '...
0558              num2str(yCoorMin + (i-1)*ydim) 'm'],...
0559             'FontName','Times New Roman',...
0560             'FontSize',14);
0561      colormap(hot)
0562      set(gca,'DataAspectRatio',[zdim xdim 1])
0563  end
0564  
0565 end

Generated on Fri 30-Jul-2010 16:56:05 by m2html © 2003