plotting cross_section using IDL

pro cross_section, input_arr, lon, lat, lonp, latp, res, $
csection = cross, xaxis = xx
;+
;NAME:
; cal_cross
;PURPOSE:
; calculate 2d cross section (radial - altitude)
; from 3d field (nx, ny, nz) and 2d lon, lat (nx, ny)
;INPUT:
; 1) two long/lat points at A and B
; 2) res ( A-X )
;NOTE
;(horizontal plain)
; A
; |\
; | X
; | \
; | X dis
; | \
; dy | X
; | \
; | X
; | \
; C_X_X_X_X_B
; dx
;
; AB(dis) : AX(res) = CB(dx) : CX (res1)

;OUTPUT
; 2 dimensional array [num, nz]
; a---------a
; | |
; | |
; | |
; A_X_X_X_X_B
; x-axis in km
;
; cross_section.pro , v 1.0 by M.-S. Park (mpark@nps.edu)
;_________________________________________________________
check_dim = size(input_arr)
if check_dim(0) eq 3 then begin
nz = check_dim(1)
nz = 31
; 1. finding lon, lat array for cross section
; (1) finding AB (dis), dx, dy in km
dlon = lonp(1) - lonp(0)
dlat = latp(1) - latp(0)
if dlon lt 0. then dlon = -dlon
if dlat lt 0. then dlat = -dlat
dx = deg2km (dlon, lonp(0), latp(0))
dy = deg2km (dlat, lonp(0), latp(0))
print, dx(0), ' km along longitude.'
print, dy(1), ' km along latitude.'
dis = (dx(0)^2. + dy(1)^2.)^(0.5)
print, dis, ' km is dis between A and B'
; (2) finding res1 in km
xx = [ 0., dx(0)]
yy = [dy(1), 0.]
res1 = res * dx(0) / dis
; (3) finding num & making cross array
num = fix(dx(0)/res1)
cross = fltarr(num, nz)
in = where(min(latp) eq latp)
; (4) res1 in km to reg1_deg.
imsi = km2deg(res1,lonp(in), latp(in))
res1_deg=imsi(0)
; (5) calculating slope
reg = linfit(lonp, latp, /double)
; (6) making lon, lat array for cross section
xx = findgen(num)*res1
lon_arr = findgen(num)*res1_deg + lonp(0)
lat_arr = fltarr(num)

for inum = 0, num - 1 do begin
lat_arr (inum) = reg(1) * lon_arr(inum) + reg(0)
endfor
; 2. finding the nearest points
ind1 = intarr(2, num) ;the first nearestpoint
ind2 = intarr(2, num) ;the second nearestpoint

for inum = 0, num - 1 do begin
pnt = Nearestpoint (lon, lat, lon_arr(inum), lat_arr(inum))
print, 'nearest index', pnt
lon_imsi = lon
lat_imsi = lat
dx = lon_arr(inum) - lon(pnt(0), pnt(1))
dy = lat_arr(inum) - lat(pnt(0), pnt(1))
; print, pnt, pnt_2
for iz = 0, nz - 1 do begin
cross(inum, iz) = input_arr (pnt(0),pnt(1),iz)
endfor
endfor
;return
endif ; some conditions for checking input
end