;+
function plot_thetas, im, xc, yc, th, HDR=hdr, INTRVAL=intrval, RANGE=range, SDIR=sdir, ROOT=root, SAVE=save, PIXELS=pixels, OCC=occ, RSUN=rsun, OVERPLOT=overplot, YTITLE=ytitle, XRANGE=xrange, YRANGE=yrange, THICKNESS=thickness
;
; Purpose:
; generates values of plots starting at a given center for different thetas
;
; INPUTS:
;  im:   	2D array
;  xc, yc:	Center of plot in IDL pixels
;  th:		(first) angle to plot
;
; KEYWORDS
;  INTRVAL:	Interval of Thetas to test (ie, 2 means plot every other angle)
;  HDR:		LASCO or EIT image header (structure)
;  RANGE:	range of angles to test (Degrees). Default is 180.
;  SDIR:		Directory to save plots in
;  ROOT:		Identifier for plot filenames
;  RSUN:	Label in units of solar radii from center (needs HDR set)
;  START:	First angle to plot (default is 0, defined as right equator)
;  OCC:		Use occulter center instead of sun center (if HDR is set)
;  RMAX:	Length of each ray computed; default is sqrt(naxis1)/2
;  OVERPLOT:	Use oplot instead of plot
;  YTITLE:	Units of values in im (string)
;
; RESTRICTIONS:
;  Must have window 0 for image and window 2 for plot open
;
; Written by N. Rich, NRL/Interferometrics
;
; Modified
;  8/21/02, N. Rich - Make general rather than using header for center
;
; 08/29/02 @(#)plot_thetas.pro	1.1
;
;-

IF keyword_set(SDIR) THEN dir=sdir ELSE dir='.'
IF keyword_set(ROOT) THEN rout = root ELSE rout = 'ptheta_'
IF keyword_set(RANGE) THEN rang = range ELSE rang=180.
IF keyword_set(INTRVAL) THEN intrv=intrval ELSE intrv=10
IF keyword_set(THICKNESS) THEN thick=thickness ELSE thick=2
sz=size(im)
maxax=(sz[1]>sz[2])
;maxim = MAX(ABS(im))
;maxim = 2500
maxim=0.5
IF keyword_set(RMAX) THEN width = rmax ELSE width = FIX(maxax*sqrt(2)/2)
help,width
n=FIX(rang/intrv)
r_all = fltarr(width,n+1)
;bias = get_bias(hdr)
;byc = interpolate(byc,x)
xax = findgen(width)
xtit = 'Pixels'

IF keyword_set(HDR) THEN BEGIN
	IF keyword_set(OCC) THEN BEGIN
   		center=occltr_cntr(hdr)
   		xc=center(0)
   		yc=center(1)
	ENDIF ELSE BEGIN
   		center = get_sun_center(hdr)
   		yc = center.ycen
   		xc = center.xcen
	ENDELSE
	arcs = get_sec_pixel(hdr)*1024/width
	divperpix = arcs/intrv
	date = hdr.date_obs
	yymmdd=strmid(date,2,2)+strmid(date,5,2)+strmid(date,8,2)
	solar_ephem,yymmdd,radius=radius,/soho
	;x = findgen(divperpix*width)
	;x = x/divperpix
	;help,x
	IF keyword_set(RSUN) THEN BEGIN
		xax = (xax-xc)*arcs/(radius*3600)
		xtit = 'Rsun'
	ENDIF
ENDIF

print,xc,yc
print,xax(0),xax(width-1)
;loadct,0
;ytit=''

;window,0,xsiz=sz[1],ysiz=sz[2]
wset,0
;tvscl,hist_equal(im)

;th = 0*intrv
titltxt = 'Point to Point profile at Theta='+string(th,format='(i3)')+' deg'
print,titltxt
;ray = im(xc-width/2:xc+width/2-1,yc)
p2 = [xc+width*cos(th*!pi/180.),yc+width*sin(th*!pi/180.)]
p1 = [xc,yc]
ray = get_profilem(im,p1(0),p1(1),p2(0),p2(1))
plots,[p1(0),p2(0)],[p1(1),p2(1)],/device
;window,2,xsiz=800,ysiz=600	;,/pixmap
wset,2
IF keyword_set(OVERPLOT) THEN oplot, xax, ray, thick=thick ELSE $
   plot,xax,ray,thick=thick, /yticklen, yticks=8, ytitle=ytitle, xtitle=xtit, title=titltxt, xrange=xrange, yrange=yrange
;stop
IF keyword_set(SAVE) THEN BEGIN
	graph = tvrd()
	gifname = dir+'/'+rout+string(th,format='(i3.3)')+'.gif'
   	write_gif,gifname,graph
ENDIF
r_all(*,0)=ray

IF keyword_set(INTRVAL) THEN BEGIN
FOR i=1,n DO BEGIN
   th = i*intrv
   titltxt = 'Point to Point profile at Theta='+string(th,format='(i3)')+' deg'
   print,titltxt
   ;p1 = [xc+500*cos(th*!pi/180.),yc+500*sin(th*!pi/180.)]
   ;p2 = [xc+500*cos((th+180)*!pi/180.),yc+500*sin((th+180)*!pi/180.)]
   p2 = [xc+width*cos(th*!pi/180.),yc+width*sin(th*!pi/180.)]
   ray = get_profilem(im,p1(0),p1(1),p2(0),p2(1))
   wset,0
   plots,[p1(0),p2(0)],[p1(1),p2(1)],/device
   testim=tvrd()
   ;sdiff = CONGRID(ray,width)
   ;window,xsiz=800,ysiz=600	;,/pixmap
   window,2
   ;tv,rebin(testim,512,512)
   wait,1
   plot,xax,ray,thick=thick, /yticklen, yticks=8, ytitle=ytitle, xtitle=xtit, title=titltxt, xrange=xrange, yrange=yrange
;stop
   graph = tvrd()
   gifname = dir+'/'+rout+string(th,format='(i3.3)')+'.gif'
   IF not(keyword_set(NOSAVE)) THEN write_gif,gifname,graph
   r_all(*,i)=sdiff
ENDFOR
return, r_all
ENDIF

;wset,3
;testim = tvrd()
;wdelete
;window,xsiz=hdr.naxis1,ysiz=hdr.naxis2
;tv,testim
return, reform(r_all[*,0])

end