* Script that takes a successive list of points, representing the 
* bounds of some region (presumably a limited-area grid) and then
* plots them as a polygon on the current figure.
* Written by Lucas Harris, Lucas.Harris@noaa.gov April 2011

* This routine is based off of basemap.gs at 
* ftp://grads.iges.org/grads/scripts/basemap.gs
* written by Jennifer M. Adams

* Information about string manipulation commands is at 
* http://www.iges.org/grads/gadoc/script.html#intrinsic

function main(args)
*function main(file, color)

file=subwrd(args,1)
color=subwrd(args,2)

if (color='')
  'set line 1 1 6'
else
  'set line 'color' 1 6'
endif

* Make sure there's a plot already drawn
'q gxinfo'
line5 = sublin(result,5)
line6 = sublin(result,6)
xaxis = subwrd(line5,3)
yaxis = subwrd(line5,6)
proj  = subwrd(line6,3)
if (xaxis = 'None' | yaxis = 'None') 
  say 'You must display a variable before using draw_grid'
  return
endif

* Get the image edges for clipping
'q gxinfo'
line3 = sublin(result,3)
line4 = sublin(result,4)
xd1 = subwrd(line3,4)
xd2 = subwrd(line3,6)
yd1 = subwrd(line4,4)
yd2 = subwrd(line4,6)
'set clip 'xd1' 'xd2' 'yd1' 'yd2

* Get the lat/lon range of the current dimension environment
  'q dims'
  line1 = sublin(result,2)
  line2 = sublin(result,3)
  minlon = subwrd(line1,6)
  maxlon = subwrd(line1,8)
  minlat = subwrd(line2,6)
  maxlat = subwrd(line2,8)

* Start reading lines. We want to read two lines from the file, 
* convert those lines (in lat-lon coordinates) into xy 
* coordinates, and then plot a line between the two.

* Read the first record from the polygon file
result = read(file)
startline = sublin(result,2)
rc = sublin(result,1)
rc = subwrd(rc,1)
if (rc!=0)
  say 'Error reading 'file
  return
endif
lon1 = subwrd(startline,1)
lat1 = subwrd(startline,2)

result = read(file)
endline = sublin(result,2)
lon2 = subwrd(endline,1)
lat2 = subwrd(endline,2)

* Convert from lat-lon to xy
'q w2xy ' startline
x1 = subwrd(result,3)
y1 = subwrd(result,6)

'q w2xy ' endline
x2 = subwrd(result,3)
y2 = subwrd(result,6)

if ((lon1-lon2) < 180 & (lon2-lon1) < 180)
 'draw line 'x1' 'y1' 'x2' 'y2
endif

* Read subsequent records
* We would also like that if a blank line is found, skip it, and 
* the next point starts a new line; i.e. there is a break in the curve
 continueline=1
 flag = 1
 while (flag)

   while(1)
     x1 = x2
     y1 = y2
     lon1 = lon2
     lat1 = lat2

     result = read(file)
     rc = sublin(result,1)
     rc = subwrd(rc,1)
     if (rc!=0)
       flag = 0
       break
     endif
     endline = sublin(result,2)
     lon2 = subwrd(endline,1)
     lat2 = subwrd(endline,2)
     if (lon='' | lat2='')
	continueline=0
	break
     endif
*     if (sublin!=0)
*       continue
*     endif

     'q w2xy ' endline
     x2 = subwrd(result,3)
     y2 = subwrd(result,6)

     if ((lon1-lon2) < 180 & (lon2-lon1) < 180 & continueline=1)
     'draw line 'x1' 'y1' 'x2' 'y2
     else
*        say lon1' 'lon2
       	continueline=1
     endif

   endwhile
 endwhile



*exit for now

*close(file)
*return