radius problem in viewshed

i m using Gtopo maps and in plotting for viewshed analysis i m getting radius problem.getting abrupt answers. > thanks for any help or reference > Hussain

Its working but there is some problem with the radius input methods.

clc figure % hold on % axesm('MapProjection','mercator','FLatlim',[23 38],'FLonlim',[62 78]); latlim = [31 33]; lonlim = [69 73]; [Z,refvec] = gtopo30('C:\GTOPO30',1,latlim,lonlim); worldmap(Z,refvec); % demcmap(Z); lighting phong; % grid off; % meshm(Z,refvec,size(Z),Z); % daspectm('m',12); % camlight; waypoints = [23.63,68.137;23.665,68.1419;23.9,68.27; 23.973 ,68.36; 23.95 ,68.755; 24.35, 68.755; 24.236, 68.8675; 24.309, 69.10667; 24.265, 69.619; 24.1519, 70.0713; 24.309, 70.13; 24.4408, 70.6122; 24.236, 70.6122; 24.2655, 70.885; 24.3677, 71.036; 24.616, 70.9922; 24.674, 71.0947; 25.3347, 70.0116; 25.7878, 70.64; 25.6561, 70.4075; 25.963, 70.07138; 26.2844, 70.20305; 26.41, 70.144; 26.6641, 69.706; 26.7958, 69.443; 27.1755, 69.603; 27.643, 70.056; 27.8475, 70.159; 28, 70.06; 27.7013, 70.7878; 27.86222, 71.3869; 27.933, 71.8986; 28.2127, 72; 28.7825, 72.351; 29.0308, 72.9508; 29.9075, 73.403; 30.2, 73.944; 30.389, 73.8277; 31.07611, 74.588; 31.458, 74.6025; 31.795, 74.5294; 32.2625, 75.348; 32.466, 74.7633; 32.773, 74.6611; 33.0366, 74.2811; 33.212, 74; 33.3288, 74.1786; 33.665, 73.959; 34, 74.2811; 34.1033, 73.9302; 34.3225, 74.01777; 34.366, 73.6816; 34.658, 73.974; 34.731, 74.5441; 34.541, 75.333; 34.4538, 75.655; 34.6147, 76.3566; 34.7313, 76.8388; 35.4911, 77.83277; 35.798, 76.59; 35.656, 76.3858; 36.0608, 75.9619; 36.4408, 75.976; 36.7333, 75.40638; 37, 75.7633; 36.7913,74.0325; 36.7475,72.364; 36.046, 71.1822; 35.827, 71.416; 35.169, 71.65; 34.526, 71.0652; 34.43916, 70.992; 34.278, 71.182; 33.951, 70.904; 34.044, 69.925; 33.898, 69.88; 33.679, 70.133; 33.489, 70.173; 33.299, 70.29; 33.022, 69.525; 32.759, 69.413; 32.4525, 69.2765; 31.926, 69.3405; 31.605, 68.91666; 31.721, 68.5219; 31.7511, 68.0836; 31.5466, 67.58666; 31.3422, 67.8058; 31.269, 67.0311; 30.947, 66.401; 30.17055, 66.285; 29.951, 66.358; 29.849, 66.241; % for boundry and additional info (this needs layer) 29.542, 65.014; 29.483, 64.1808; 29.381, 64.0933; 29.4252, 62.485; 29.8783, 60.863; 29.381, 61.316; 29.0016, 61.535; 28.636, 61.886; 28.475, 62.2516; 28.358, 62.5; 28.256, 62.587; 28.256, 62.806; 27.979, 62.7486; 27.4969, 62.965; 27.219, 62.777; 27.263, 63.215; 27.1025, 63.303; 26.839, 63.289; 26.839, 63.1836; 26.6058, 63.1725; 26.5911, 62.7194; 26.518, 62.339; 26.284, 62.226; 26.284, 62.076; 26.226, 61.871; 25.99, 61.827; 25.743, 61.856; 25.67, 61.969; 25.145, 61.5938;]; [lttrk,lntrk] = track('rh',waypoints,'degrees');

comet3m(lttrk,lntrk,8000); plot3m(lttrk,lntrk,'k','Erasemode','none','LineWidth',2);

lat1=32; lon1=71; % actualradius = almanac('earth','radius','m'); actualradius=637100; apparentradius =actualradius; oalt=10; talt=300; R=refvec2mat(refvec,size(Z)); % Re=8000; %for one lat long [visgrid,visleg] = viewshed(Z,refvec,lat1,lon1,oalt,talt,'AGL','AGL',actualradius,apparentradius); % axesm('globe','geoid',almanac('earth','wgs84','meters')) meshm(visgrid,visleg,size(Z),Z); caxis auto; colormap(cool(2)); lighting phong; daspectm('m',5); camlight; h = lcolorbar({'obscured','visible'}); set(h,'Position',[.875 .45 .02 .1])