for ir = 1:21
    for ith = 1:61
        r = 0.5+(ir-1)*0.025;
        th = (ith-1)*(pi/30);
        x(ir,ith) = r*cos(th);
        y(ir,ith) = r*sin(th);
        c(ir,ith) = r*cos(3*th);
    end
end
pcolor(x,y,c)
shading interp
hold on
contour(x,y,c,'k-')
colorbar
axis equal tight