; NAME :	dispr
; PURPOSE :
;	plot a sphere and indicate the observed region on it
; CATEGORY :
; 	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	dispr, mparam
; INPUTS :
;	mparam
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	none
; OUTPUTS :
;	plot a sphere and a rectangle indicating the observed region
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	X-window environment is assumed
; PROCEDURE :
;
; MODIFICATION HISTORY :
;	T.Sakurai, 1992
;
;------------------------------------------------------------------------
;
function matrix_r, mparam
;
	prad = mparam.prad
	b0rad = mparam.b0rad

	armat = fltarr(3,3)
	armat(0,0)= cos(prad)
	armat(0,1)=-sin(prad)*cos(b0rad)
	armat(0,2)= sin(prad)*sin(b0rad)
	armat(1,0)= sin(prad)
	armat(1,1)= cos(prad)*cos(b0rad)
	armat(1,2)=-cos(prad)*sin(b0rad)
	armat(2,0)= 0.0
	armat(2,1)= sin(b0rad)
	armat(2,2)= cos(b0rad)
	return, armat
end

pro dispr, mparam
;
	xp = fltarr(1000)
	yp = fltarr(1000)
	pi = 3.141592
	armat = matrix_r(mparam)

;	wdelete, 0, 0 	; ? wdelete,0 gives no action

	wsize = 350
      window, 2, title='IDL 2 : MAGPACK2 display', $
		xsize=wsize, ysize=wsize
      setcolor, bandwidth, w_start, rb_start, g_start, color8

            xcorner = 0
            ycorner = 0
; xcorner,ycorner = lower left corner of plot window in device coordinates

; lower left and upper right corners of contour plot in device coordinates
            mapposition = [xcorner, ycorner, $
                        xcorner+wsize -1, ycorner+wsize -1]

	rsp = mparam.rsp
	rspmargin = 0.95*rsp
	plot, [-rspmargin,rspmargin],[-rspmargin,rspmargin],/nodata, $
		position = mapposition, /device, $
		xstyle=4, ystyle=4

	xyouts, -rsp,rsp+5,'!6'+mparam.title

	eps=-1.0e-2
	nc=200
	xp = rsp*cos(2.0*pi*indgen(nc)/(nc-1)) 
	yp = rsp*sin(2.0*pi*indgen(nc)/(nc-1))
	oplot,xp,yp, color = color8.yellow
;
	for i=1,17 do begin
		th=10.0/180.0*pi*i
		np=fix(nc*sin(th))
		if(np lt 6) then np=6
		ph = 2.0*pi*indgen(np)/(np-1) - 0.5*pi
		ux1=rsp*sin(th)*cos(ph)
		ux2=rsp*cos(th)
		ux3=rsp*sin(th)*sin(ph)
		xp=armat(0,0)*ux1+armat(0,1)*ux2+armat(0,2)*ux3 
		yp=armat(1,0)*ux1+armat(1,1)*ux2+armat(1,2)*ux3 
		zp=armat(2,0)*ux1+armat(2,1)*ux2+armat(2,2)*ux3

		visible = where( zp ge eps, count )
		if count gt 0 then oplot, xp(visible), yp(visible), $
			color = color8.yellow
	endfor
;
	for i=1,36 do begin
		ph=10.0/180.0*pi*i
		dth=1.0
		if( (i mod 9) eq 0) then dth=0.0
		wth=18.0-2.0*dth
		nt=nc/2
		th = pi*(wth*indgen(nt)/(nt-1) +dth)/18.0
		ux1=rsp*sin(th)*cos(ph)
		ux2=rsp*cos(th)
		ux3=rsp*sin(th)*sin(ph)
		xp=armat(0,0)*ux1+armat(0,1)*ux2+armat(0,2)*ux3 
		yp=armat(1,0)*ux1+armat(1,1)*ux2+armat(1,2)*ux3 
		zp=armat(2,0)*ux1+armat(2,1)*ux2+armat(2,2)*ux3

		visible = where( zp ge eps, count )
		if count gt 0 then oplot, xp(visible), yp(visible), $
			color = color8.yellow
	endfor
;
	xs = mparam.xs
	ys = mparam.ys
	xe = xs + mparam.dx*(mparam.nx -1)
	ye = ys + mparam.dy*(mparam.ny -1)
	oplot,[xs,xe,xe,xs,xs],[ys,ys,ye,ye,ys]

	wset, 0
	return
end