%
%   BARE1G_RZ.M   Plot Spatial Flux and Current Distribution 
%                   for a 1-G RZ Bare Cylindrical Reactor 
%
%   This demo displays the flux and current profiles for a bare homogeneous critical
%   finite cylindrical reactor.  For convenience, the extrapolation distance is 
%   assumed to be small relative to the reactor geometry.  In addition, the 
%   fundamental mode flux is normalized to a maximum value of unity and the 
%   diffusion coeff in the current equation is also unity. 
%
%   Code written by J. R. White, UMass-Lowell (Jan. 2015)
%

      clear all, close all, nfig = 0;
%
%   calc normalized fluxes & currents 
%    --> let R = 1 m, then H = 2*Dia = 4R = 4 m
%    --> let A = 1 in the flux equation  and  let DA = 1 in the current equation 
      R = 1;  H = 4*R;  a = 2.4048;  
      npts = 21;  r = linspace(0,R,npts);   z = linspace(0,H/2,npts);
      [rr,zz] = meshgrid(r,z);
      Xr = besselj(0,a*rr/R);  Yz = cos(pi*zz/H);        % radial & axial components
      flxrz = Xr.*Yz;                                    % flux profile (A = 1)
      currzr = (a/R)*besselj(1,a*rr/R).*cos(pi*zz/H);    % r-directed current
      currzz = (pi/H)*besselj(0,a*rr/R).*sin(pi*zz/H);   % z-directed current
      mcurrz = sqrt(currzr.*currzr + currzz.*currzz);    % current magnitude (DA = 1)
%
%   plot fluxes
      nfig = nfig+1;     figure(nfig);
      plot(r,Xr(1,:),'r-','LineWidth',2),grid on
      title('Normalized Radial Profile in Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Normalized Flux')
%
      nfig = nfig+1;     figure(nfig);
      plot(z,Yz(:,1),'b--','LineWidth',2),grid on
      title('Normalized Axial Profile in Bare RZ Reactor Model')
      xlabel('Axial Distance (m)'),ylabel('Normalized Flux')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,flxrz),view(45,30),shading interp, colorbar, axis image
      title('Fundamental Mode Flux for Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized Flux')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,flxrz),view(2),shading interp, colorbar, axis image
      title('Fundamental Mode Flux for Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized Flux')
%
%   plot currents
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,currzr),view(45,30),shading interp, colorbar, axis image
      title('r-Directed Current in Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized r-Directed Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,currzr),view(2),shading interp, colorbar, axis image
      title('r-Directed Current in Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized r-Directed Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,currzz),view(45,30),shading interp, colorbar, axis image
      title('z-Directed Current Distribution for RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized z-Directed Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,currzz),view(2),shading interp, colorbar, axis image
      title('z-Directed Current Distribution for RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized z-Directed Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,mcurrz),view(45,30),shading interp, colorbar, axis image
      title('Current Magnitude in Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      surf(rr,zz,mcurrz),view(2),shading interp, colorbar, axis image
      title('Current Magnitude in Bare RZ Reactor Model')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
      zlabel('Normalized Current')
%
      nfig = nfig+1;     h = figure(nfig);  set(h,'Renderer','Zbuffer')
      vc = [.3 .5 .7 .8 1.0 1.2 1.3 1.38];
      [c,h] = contour(rr,zz,mcurrz,vc); clabel(c); hold on
      set(h,'Linewidth',2)
      quiver(rr,zz,currzr,currzz), hold off, axis image
      title('Current Contour Lines & Directions (RZ Model)')
      xlabel('Radial Distance (m)'),ylabel('Axial Distance (m)')
%
%   end of simulation
