pro ps_profiles, i_obs,q_obs,u_obs,v_obs,have_un,i_mag,i_fit,q_fit,u_fit, $
	v_fit,nounnofit,xs,ys,encapsulated=encapsulated_, noImag=noImag_
;------------------------------------------------------------------------------
;
;	procedure:  ps_profiles
;
;	purpose:  plot observed, fitted Stokes profiles from AZAM for one 
;			  raster point.  Plot to PostScript files.
;
;	inputs:	iobs,qobs,uobs,vobs - observed stokes profiles
;		imag - i_magnetic fitted profile
;		ifit,qfit,ufit,vfit - fit stokes profiles
;		have_un	- flag, if non-zero then unno fit profiles at this pt
;		xs,ys - coordinates in (reduced) image

;
;------------------------------------------------------------------------------

				    ;Check for encapsulated PostScript.
encapsulated = keyword_set( encapsulated_ )
extd = (['ps','eps'])(encapsulated)
				    ;Check if i-magnetic plot suppressed.
noImag = keyword_set( noImag_ )

;  common area for usage in azam_spectra.pro
common com_azamspectra, xllo,xuro,yllo,yuro,full_directory, $
	iwstrt,iwend,dispers,wavlngth,basewav,etaratio,vbroad, $
	axsz,ntran

;  calculate the pixel ranges for wavelength cutouts
nwrng = n_elements(iwstrt)
;  nlinr = ranges of pixels for output array for each spectrum line
nlinr = intarr(nwrng,2)
;  nwave = total number of wavelengths to output
nwave = 0
for kk = 0,nwrng-1 do begin
	nwlin = (iwend(kk)-iwstrt(kk) + 1)
        nlinr(kk,0) = nwave
        nlinr(kk,1) = nlinr(kk,0) + nwlin - 1
        nwave = nwave + nwlin
endfor

;  set dispersion in angstroms
disper = dispers*.001

;  calculate the wavelengths for the observed range based on
;  the nominal wavelengths
nwobs = n_elements(i_obs)	;  number of observed wavelength pixels
wobs = dblarr(nwobs)
;  starting value of wavelength for observations
wobstrt = basewav(0) - disper(0)*float(iwstrt(0))
;  below assumes that a continuous wavelength range is observed
wobs = wobstrt + disper(0)*findgen(nwobs)
;  correct the wavelengths for the remaining lines
for kk = 1,nwrng-1 do begin
	for ii = iwstrt(kk),iwend(kk) do begin
		wobs(ii) = basewav(kk) + disper(kk)*float(ii-iwstrt(kk))
	endfor
endfor
;  calculate the nominal wavelengths for the unno profiles
wunno = dblarr(nwave)
for kk = 0,nwrng-1 do begin
	for ii = nlinr(kk,0),nlinr(kk,1) do begin
		wunno(ii) = basewav(kk) + disper(kk)*(float(ii)-nlinr(kk,0))
	endfor
endfor

;  get the original data pixel coordinates
xp = xs + xllo & yp = ys + yllo


;Save system stuff.
sav_p=!p  &  sav_x=!x  &  sav_y=!y  &  sav_z=!z  &  sav_o=!order
tvlct, sav_r, sav_g, sav_b, /get
sav_w=!d.window  &  sav_n=!d.name



;  Begin PostScript stuff
;  set character, line thickness
thick = 5.
!x.thick=thick
!y.thick=thick
cthick = 3.
charthick = cthick
  ;
  ;       Set the device type to high-resolution PostScript.
  ;
  set_plot, 'ps'
;  set for portrait mode
device, bits_per_pixel=8, /inches,scale_factor=1.0, $
	;file='plotqvprofs.'+extd, scale_factor=1.0, /inches, $
	file=strcompress(/remove_all,'profiles_'+string(xp)+'_'+string(yp)+'.'+extd), $
	xoffset=0.0, yoffset=0.0, xsize=8.5, ysize=11.0, /portrait, color=1, $
	/times, /bold, encapsul=encapsulated


;  set the wavelength annotations and ticks
;  get the central wavelength of the observed range to the nearest 0.5 A
wavmid = double(round(2.*mean(wobs,/double))/2.)
print,'mid wavelength: ',wavmid
; find ticks at lower wavelengths
tkmin = wavmid & tkmax = wavmid
ntmin = 0
for ii = 0,100 do begin
	tkmin = tkmin - 0.5
	ntmin = ntmin+1
	if tkmin lt min(wobs) then begin
		tkmin = tkmin + 0.5
		ntmin = ntmin-1
		goto,jump55
	endif
endfor
jump55: for ii = 0,100 do begin
	tkmax = tkmax + 0.5
	if tkmax gt max(wobs) then begin
		tkmax = tkmax - 0.5
		goto,jump56
	endif
endfor
jump56:
;  set wavelength tick marks
nxt = round(double(tkmax-tkmin)/0.499)
tickx =  wavmid - ntmin*0.5 + findgen(nxt+1)*0.5
;  get labels for ticks
xtn = strmid(strcompress(string(tickx-wavmid),/remove_all),0,4)
;  set middle of wavelength to full value in A
;  This presumes the wavelength is 4 characters
xtn(ntmin) = strmid(strcompress(string(wavmid),/remove_all),0,6)

;Load black & white color table.
loadct, 0
ncidx = 254
;  set special colors
black  = 0        &  tvlct,   0,   0,   0, black
white  = ncidx-1  &  tvlct, 255, 255, 255, white
yellow = ncidx-3  &  tvlct, 255, 255,   0, yellow
red    = ncidx-4  &  tvlct, 255,   0,   0, red
green  = ncidx-5  &  tvlct,   0, 255,   0, green
blue   = ncidx-6  &  tvlct,   0,   0, 255, blue

; some middle index
sz = 300

;;Set plot charsize.
charsize = 1.0

;  specify wavelength range
rangew = [wobs(0),wobs(nwobs-1)]

				    ;Factor to allow margin.
	fm = 1.05
				    ;Plot unsolved point.
	if not have_un then begin
		;		    ;Do plotting in pixmap window.
		;wset, wp
				    ;Plot i profile.
		!p.region = [ 0.0, 0.5, 0.5, 1.0 ]
		max_i = max(i_obs)
		mx = fm*max_i
		;plot, [0,numx-1], [0,mx], ystyle=1, /nodata $
		plot, rangew, [0,mx], ystyle=1, /nodata $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
				
		plots, wobs, i_obs, color=black

				    ;Plot q profile.
		!p.region = [ 0.5, 0.5, 1.0, 1.0 ]
		max_q = max(abs(q_obs)) > max(abs(u_obs)) 
		;if fix_q ne 0. then  mx=fix_q  else  mx=max_q
		mx=max_q
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		plots, wobs, q_obs, color=black

				    ;Plot u profile.
		!p.region = [ 0.0, 0.0, 0.5, 0.5 ]
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		plots, wobs, u_obs, color=black

				    ;Plot v profile.
		!p.region = [ 0.5, 0.0, 1.0, 0.5 ]
		max_v = max(abs(v_obs))
		;if fix_v ne 0. then  mx=fix_v  else  mx=max_v
		mx=max_v
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		plots, wobs, v_obs, color=black

				    ;Print some information.
;  print out coordinates from start of full map
		xyouts, 0.40, 0.49, color=black, /normal, 'Image Pixel '+ $
		'['+stringit(xllo+xs)+','+stringit(yllo+ys)+']', charsize=1.3

	end




				    ;Plot profiles & Unno fit.
	if have_un then begin
				    ;Plot i profiles with Unno fit.
		!p.region = [ 0.0, 0.5, 0.5, 1.0 ]
		max_i = max(i_obs) > max(i_fit)
		mx=max_i
		;plot, [0,numx-1], fm*[0,mx], ystyle=1, /nodata $
		plot, rangew, fm*[0,mx], ystyle=1, /nodata $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		for iline = 0,nwrng-1 do begin
;  plot fitted i_mag in blue
			if(noImag) then goto, jump66
			plots, wunno(nlinr(iline,0):nlinr(iline,1)), $
				i_mag(nlinr(iline,0):nlinr(iline,1)), color=blue
			jump66:
			plots, wunno(nlinr(iline,0):nlinr(iline,1)), $
				i_fit(nlinr(iline,0):nlinr(iline,1)), color=red
		endfor
;  bypass magnetic only I-profile plotting
		plots, wobs, i_obs, color=black

				    ;Plot q profile with Unno fit.
		!p.region = [ 0.5, 0.5, 1.0, 1.0 ]
		max_q = max(abs(q_obs)) > max(abs(u_obs)) $
	 	      > max(abs(q_fit)) > max(abs(u_fit)) 
;		if twoComp then $
;		max_q = max_q $
;	 	      > max(abs(q_1st)) > max(abs(u_1st)) $
;	 	      > max(abs(q_2nd)) > max(abs(u_2nd))
		;if fix_q ne 0. then  mx=fix_q  else  mx=max_q
		mx=max_q
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		for iline = 0,nwrng-1 do begin
			plots, wunno(nlinr(iline,0):nlinr(iline,1)), $
				q_fit(nlinr(iline,0):nlinr(iline,1)), color=red
		endfor
		plots, wobs, q_obs, color=black

				    ;Plot u profile with Unno fit.
		!p.region = [ 0.0, 0.0, 0.5, 0.5 ]
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		for iline = 0,nwrng-1 do begin
			plots, wunno(nlinr(iline,0):nlinr(iline,1)), $
				u_fit(nlinr(iline,0):nlinr(iline,1)), color=red
		endfor
		plots, wobs, u_obs, color=black

				    ;Plot v profile with Unno fit.
		!p.region = [ 0.5, 0.0, 1.0, 0.5 ]
		max_v = max(abs(v_obs)) > max(abs(v_fit))
;		if twoComp then $
;		max_v = max_v > max(abs(v_1st)) > max(abs(v_2nd))
		;if fix_v ne 0. then  mx=fix_v  else  mx=max_v
		mx=max_v
		plot, rangew, fm*[-mx,mx], ystyle=1, /nodata, /noerase $
		, charsize=charsize, background=white, color=black, $
		xminor=5, xtickv = tickx, xtickname=xtn, xticks=nxt
		for iline = 0,nwrng-1 do begin
			plots, wunno(nlinr(iline,0):nlinr(iline,1)), $
				v_fit(nlinr(iline,0):nlinr(iline,1)), color=red
		endfor
		plots, wobs, v_obs, color=black

				    ;Print some information.
;  print out coordinates from start of full map
		xyouts, 0.40, 0.49, color=black, /normal, 'Image Pixel '+ $
		'['+stringit(xllo+xs)+','+stringit(yllo+ys)+']', charsize=1.3

;  endif have_un
	end



device, /close_file


set_plot, 'x'

;Restore system stuff.
set_plot, sav_n
if sav_n eq 'X' and sav_w gt -1 then  wset,sav_w
tvlct, sav_r, sav_g, sav_b
!p=sav_p  &  !x=sav_x  &  !y=sav_y  &  !z=sav_z  &  !order=sav_o




return

end