;+
; NAME:
;	SMEI_SEQUENCE::MAKE_PROFILE
;
;
; PURPOSE:
;	Generate a radial or circumferential profile set for a SMEI
;	sequence.
;
;
; CATEGORY:
;	SMEI SEQUENCE
;
;
; CALLING SEQUENCE:
;	profile = seqref->make_profile(npa, nel, pamin, pamax)
;		or
;	profile = seqref->make_profile(npa, nel, elmin, elmax, $
;			/circumferential)
;
;
; INPUTS:
;	nel	int	The number of bins in the elongation axis.
;	npa	int	The number of bins in the position angle axis.
;	pamin	float	The lowest position angle to use
;	pamax	float	The highest position angle to use
;			 In future if pamin >
;			pamax, then the sector will be assumed to
;			cross north, but for the moment pamin must be
;			less than pamax.
;	elmin	float	The lowest elongation angle to use.
;	elmax	float	The highest elongation angle to use.
;
;
; KEYWORD PARAMETERS:
;	/circumferential	If set, then make a ring all around
;				the image at a fixed elongation range.
;	pa	float	output	A named variable to hold the position
;				angles of the bin centres.
;	elong	float	output	A named variable to hold the
;				elongations of the bin centres.
;	date	double	output	A named variable to receive the dates
;				of the images.
;	/verbose		If set, then show the region during
;				the generation of the first image.
;	plane	int	input	Specify a raw image plane to use
;				(default=processed) 
;	mask	uint	input	A mask of to apply to the flag plane
;				in version 2 data. (e.g. '800'xu to
;				exclude Sun contaminated data).
;	/miss_nan		If set, then a cell with no values is
;				NaN rather than zero.
;
; OUTPUTS:
; 	profile	float	An nel x npa x nimages array holding the
; 			profiles.
;
;
;
; MODIFICATION HISTORY:
;	Original: 12/12/03; SJT
;	Add plane keyword: 12/1/04; SJT
;	Replace old clumsy REVERSE_INDEX method with pointers owing to
;	off by one per step error: 20/7/04; SJT
;	Fix memory leak: 6/7/09; SJT
;	Add mask and miss_nan keywords: 8/7/09; SJT
;-

function smei_sequence::make_profile, npa, nel, pmin, pmax, $
                      circumferential = circumferential, pa = pa, $
                      elong = elong, date = date, verbose = verbose, $
                      plane = plane, mask = mask, miss_nan = miss_nan

img = self -> get_first()
img -> mk_pa_map, elon, posa

nimg = self -> get_count()
proj = img -> get_projection()
if (proj lt 0 or proj gt 2) then begin
    smei_msg, /alert, 'Cannot create a profile for an unknown projection' 
    return, 0
endif

if keyword_set(circumferential) then begin
    ipa = floor(npa*posa/360.)
    iel = floor(nel*(elon-pmin)/(pmax-pmin))
    pa = (findgen(npa)+0.5)*360./npa 
    elong = (findgen(nel)+0.5)*(pmax-pmin)/nel + pmin
endif else begin
    if proj eq 1 then emax = 135. $
    else emax = 180.

    posa -= pmin
    locs = where(posa lt 0., nn)
    if nn ne 0 then posa[locs] += 360.
    locs = where(posa ge 360., nw)
    if nw ne 0 then posa[locs] -= 360.

    ipa = floor(npa*posa/(pmax-pmin))
    iel = floor(nel*elon/emax)
    pa = (findgen(npa)+0.5)*(pmax-pmin)/npa + pmin
    elong = (findgen(nel)+0.5)*emax/nel
endelse

if n_elements(plane) eq 0 then plane = self -> get_plane() $
else if plane lt 0 or plane gt self -> get_n_planes() then begin
    smei_msg, /alert, "Plane index out of range"
    return, 0
endif

date = dblarr(2, nimg)


if keyword_set(miss_nan) then profile = replicate(!values.f_nan, npa, $
                                                  nel, nimg) $
else profile = fltarr(npa, nel, nimg)
cells = ptrarr(npa, nel)
ncells = lonarr(npa, nel)

iimg = 0

repeat begin
    if plane eq 0 then ima = img -> get_image() $
    else ima = img -> get_raw_image(plane = plane)
    date[*, iimg] = img -> get_time()
    flags = img -> get_flag_data()
    isflag = (n_elements(flags) gt 1)

    if keyword_set(verbose) and iimg eq 0 then begin
        sz = size(ima)
        jimg = bytarr(sz[1], sz[2])
    endif
    for ii = 0l, nel-1 do begin
        for jj = 0l, npa-1 do begin
            if iimg eq 0 then begin
                cells[jj, ii] = ptr_new(where(ipa eq jj and iel eq ii, nc))
                ncells[jj, ii] = nc
            end
            if ncells[jj, ii] ne 0 then begin
                if isflag && keyword_set(mask) then begin
                    cls = ima[*cells[jj, ii]]
                    flg = flags[*cells[jj, ii]]
                    locs = where((flg and mask) eq 0, nvc)
                    if nvc eq 0 then continue
                    profile[jj, ii, iimg] = $
                      total(cls[locs])/float(nvc)
                endif else profile[jj, ii, iimg] = $
                  total(ima[*cells[jj, ii]])/ncells[jj, ii]
                if keyword_set(verbose) and iimg eq 0 then begin
                    jimg[*cells[jj, ii]] = 255b
                    tv, jimg
                endif
            endif
        endfor
    endfor
    img = img -> get_next()
    iimg++
endrep until not obj_valid(img)

ptr_free, cells

return, profile

end