
lat=-90:2.5:90;
lon=0:2.5:360;


height=3000;
olon=90; %origin longitude
olat=40;
wlon=12.5;
wlat=12.5;
rlon=0;
rlat=20;

NN=length(lat);
N=length(lon);

for n=1:N
    for nn=1:NN
        dx=min(abs([lon(n)-(olon-rlon) lon(n)-(olon+rlon)]));
        if dx <= rlon;
            dx=0;
        end

        %dx=dx*cos((pi/180)*lat(nn));
        
        dy=min(abs([lat(nn)-(olat-rlat) lat(nn)-(olat+rlat)]));
        if dy <= rlat;
            dy=0;
        end
        
        TOP(nn,n)=height*(exp(-(dx/wlon)^2-(dy/wlat)^2));
        
        
        DX(nn,n)=dx;
        DY(nn,n)=dy;
        
        clear dx dy;
    end
end

[LON, LAT]=meshgrid(lon,lat);

figure(1); clf;
pcolor(LON,LAT, TOP); colorbar;
shading interp

%figure(2); clf;
%pcolor(LON,LAT, DX); colorbar;

%figure(3); clf;
%pcolor(LON,LAT, DY); colorbar;


%%%%%%%%%%%%%%%%%%%%%%%%% now interpolate to t42;

trunc=42;

[xx,yy, field_out,spec] = spharm_trans_filt(lon, lat,TOP, trunc, 1:trunc);

sphcolor2(field_out, yy*180/pi, xx*180/pi, 12);
caxis([-80 80]); colorbar;
title({['Gauss top with wlon= ', num2str(wlon), ' and wlat = ', num2str(wlat)]; [' and abrupt spectral truncation at T',num2str(trunc), ' in meters']} , 'FontSize', 18)






        





