Nx=124;
Ny=Nx;
      Heating = 2.0;

      Rbuoy=0.02;
      Rheat=0.15;
      a1=(Rbuoy-2.0*Rheat)/(4.0*Rbuoy*Rheat);
      a2=1.0/(8.0*Rbuoy*Rheat);
      delX=0.01;
      iC=Nx/2;
      jC=Ny/2;
          for J=1:Ny
          for I=1:Nx
            iD = I-iC;
            jD = J-jC;
            rD = sqrt(iD*iD+jD*jD)*delX;
            if ( rD < (Rheat-Rbuoy) ) 
              Heat(I,J) = Heating;
            elseif (rD < (Rheat+Rbuoy)) 
              Rcool = rD - Rheat - Rbuoy;
              Heat(I,J) = Heating*(a2*Rcool^(2.0)+a1*Rcool);
            else
              Heat(I,J) = 0;
            end
          end
        end
figure
pcolor(Heat(:,:));axis square;colorbar;shading flat
figure
plot(Heat(:,Nx/2))