;+
; NAME: find_nl
; PURPOSE: This program loads a MDI magnetogram, takes a subregion for each active region, then
;            finds the neutral line of each active region, and plots it overlaid on the region.
; CATEGORY: NL identification and plotting
; CALLING SEQUENCE:      load_mag,file,&&&&&&&
; INPUTS:      file              MDI magnetogram to load
;               &&&&&&&&&&&      lower-left and upper-right corners of the subfield of view
; OUTPUTS:     a ps file
; KEYWORD:              &&&&&&
; SIDE EFFECTS: none
; RESTRICTIONS:         Need solarsoft to read in fits files, &&&&&&&&.
;                       Runs into problems when given a subfield containing
;			solar limb.  Corrections to ARsel prevent this, if used
; PROCEDURE:     This program
;                       1) Reads in the MDI magnetogram  
;                       2) Smoothes and alters as per settings
;                       3) Caluclates a potential field map assuming subfield at disk center
;                       4) Calculates smoothed positive and negative field maps of the subfield of view 
;                       5) Identifies all vertices along the 0G contour
;                       6) Suppresses weak field segments using combination of pBt, and smoothed positive/negative fields
;			7) Connects all valid NL points, stores NL information
;                       8) Overlays NL map ontop of subfield of view of MDI
; MODIFICATION HISTORY:  July 2006      David Falconer          Extracted appropiate codes, for use by Hoeksema
;				 	Natty Bokenkamp		Edited extensively, added in line connecting functions
;-								Included AR selection from NOAA database
;								Added reiterative smoothing
@alex.pro
@isth.pro
@edpt.pro
@trimcoords.pro
@cropcont.pro
pro find_nl_n, header, mdi, mdi_src, no_line, nl_xy, nl_info, edpt, coords, nl_xy_edit, nl_info_edit, nl1, pbt, grad_width=grad_width, detail=detail, thresh=thresh, pbt_thresh=pbt_thresh, minsize=minsize, sep=sep, isth_cut=isth_cut, comp=comp, nl_xy_comp=nl_xy_comp, nl_info_comp=nl_info_comp, z=z, loud=loud
 
; Get value of NAN pixels from header
	if n_tags(header) gt 0 then blank=fxpar(header.info,'BLANK') $
	else blank=fxpar(header,'BLANK')

; Get unaltered copy of image from mdi_src
	no_line=0
	mdi=mdi_src
; Set NAN pixels to 0
	mdi=mdi*(mdi ne blank)

; Smooth mdi subarray at level specified by detail keyword
	if n_elements(detail) ne 0 then begin
	if detail gt 2 then mdi=smooth(mdi, detail) else if detail eq 2 then mdi=smooth(mdi,3)
	endif

; Sets gradient width, used in creating pos and neg gradient fields
	if n_elements(grad_width) eq 0 then grad_width=10
	if grad_width le 2 then grad_width=3

; Sets z value, height at which potential field is calculated
	if n_elements(z) eq 0 then z=0

; Determine potential fields
; ALEX is code that calculates a constant alpha field
; Used here with an alpha of 0 to get a potential field
; Lati and longi will tell ALEX where the subfield is located
; Border=30 is used to give a large buffer of 0G before replicating
; See accompaning ALEX code for details.
	ALEX,8,8,0,z,mdi,BXP,BYP,BZP,0,LATI=LATI,LONGI=LONGI,BORDER=30
	pbt=(bxp^2+byp^2)^0.5   ; magnitude of simulated vector magnetogram
	pos=smooth(mdi>0,grad_width)     ; smoothed positive field
	neg=smooth(mdi<0,grad_width)     ; smoothed negative field

; Identify vertices along 0G contour.
	contour, mdi, level=[0.0], path_info=info, path_xy=coords, /path_data_coords

; Creates structure for each segment along 0G contour.
	nl_str=create_struct($
		'x',-1.,'y',-1.,$		; x and y values
                'l',-1.,'id',-1,'idx',-1L,$     ; length, curve id #, xy index
                'pBt',-1.,'pos',0.,'neg',0.)    
 				; potential field, +, and - values at midpoint
	nl_int=replicate(nl_str,n_elements(coords(0,*)))

; Load x, y, id, and idx from info and coord vars
	coordi=transpose(coords)
	for i=0, (n_elements(info) - 1) do begin
        nl_int(info(i).offset:info(i).offset+info(i).n-1).x = coordi(info(i).offset:info(i).offset+info(i).n-1,0)
        nl_int(info(i).offset:info(i).offset+info(i).n-1).y = coordi(info(i).offset:info(i).offset+info(i).n-1,1)
	nl_int(info(i).offset:info(i).offset+info(i).n-1).id = i
        for j=info(i).offset, info(i).offset+info(i).n-1 do nl_int(j).idx=j
	endfor

; Determine strength of smoothed positive, smoothed negative and potential field
; at points
	nl_int.pos=interpolate(pos,nl_int.x,nl_int.y)
	nl_int.neg=interpolate(neg,nl_int.x,nl_int.y)
	nl_int.pbt=interpolate(pbt,nl_int.x,nl_int.y)

; Sets threshholds for gradient and for magnitude of vector magnetic field
; Can be set in call
	if n_elements(thresh) eq 0 then thresh=8.
	if n_elements(pbt_thresh) eq 0 then pbt_thresh=200

; Determine good neutral lines segments
        if n_elements(where((nl_int.pbt gt pbt_thresh) and (nl_int.pos gt thresh) and (nl_int.neg lt -thresh))) gt 1 then begin
	nl1=nl_int(where((nl_int.pbt gt pbt_thresh) and (nl_int.pos gt thresh) and (nl_int.neg lt -thresh)))
	endif else if where((nl_int.pbt gt pbt_thresh) and (nl_int.pos gt thresh) and (nl_int.neg lt -thresh)) ne -1 then nl1=nl_int(where((nl_int.pbt gt pbt_thresh) and (nl_int.pos gt thresh) and (nl_int.neg lt -thresh))) else no_line=1

	if no_line eq 0 then begin

; Sets the minimum contour length in data points to be considered (call var?)
	if n_elements(minsize) eq 0 then minsize=5

; Separation variable: if a good segment is separated by this much from other
; good segments, it will be discounted
; If a group of good segments is separated by this much, they will be set
; into a separate segment
	if n_elements(sep) eq 0 then sep=30

; Selects the contours which have good neutral line segments and are longer
; than minsize, puts all nl1 points for contour i into edpt(i) structure
	lst=intarr(n_elements(info))
	edpt=replicate({pt:lonarr(n_elements(nl1)+1), sep:0, start:0, fin:0}, n_elements(info))

	for i=0, n_elements(nl1)-1 do BEGIN 
         if info(nl1(i).id).n gt minsize then lst(nl1(i).id)=lst(nl1(i).id)+1

	ENDFOR
	for i=0, (n_elements(info)-1) do if lst(i) gt 1 then begin
	if (info(i).n/3) lt sep then edpt(i).sep=info(i).n/3 else edpt(i).sep=sep
	edpt(i).start=info(i).offset & edpt(i).fin=info(i).offset+info(i).n-1
	k=0 & for j=0, (n_elements(nl1)-1) do if (nl1(j).id eq i) then begin
	edpt(i).pt(k)=nl1(j).idx & k=k+1
	endif & endif
	qp=where(lst gt 1) 
	if n_elements(qp) gt 1 then begin
	cont_info=info(qp) & edpt=edpt(qp)
	endif else if qp(0) ne -1 then begin
	cont_info=info(qp) & edpt=edpt(qp)
	endif else begin
	no_line=1 & goto, no_lines
	endelse

; Condense each edpt structure
	for i=0, n_elements(edpt)-1 do begin
	tmp=edpt(i) & edpt, tmp & edpt(i)=tmp
	endfor

; Counts the number of contours that will be needed for nl_info using endpoints
	nc=0 & for i=0, n_elements(edpt)-1 do begin
	q=where(edpt(i).pt)
	if q(0) ne -1 then begin
	if edpt(i).pt(0) ne 0 then l=n_elements(where(edpt(i).pt)) else l=n_elements(where(edpt(i).pt))+1
	if edpt(i).pt(max(where(edpt(i).pt))) eq edpt(i).fin and ((l mod 2) eq 1) then l=l-1
	if (l gt 2) then nc=nc+(l/2) else nc=nc+1
	if (edpt(i).pt(max(where(edpt(i).pt))) eq edpt(i).fin and edpt(i).pt(0) eq edpt(i).start) then nc=nc+1
	endif
	endfor

; Declares nl_info structure which will hold contour info for output
	nl_info=replicate({offset:0L,n:0L,typ:0},nc)

; Comparative NL selection is off by default
	if n_elements(comp) eq 0 then comp=0

; Comparative smoothing routine: takes information from nl_xy_comp
; and nl_info_comp variables and uses those to select the good neutral lines
; from among the potential contours and to set endpoints
	if (comp) then begin
	q=where((nl_info_comp.typ ne 2) and (nl_info_comp.n gt (max(nl_info_comp.n)/2 < 15)))

	if q(0) eq -1 then begin & no_line=1 & goto, no_lines & endif
	nl_info_comp=nl_info_comp(where((nl_info_comp.typ ne 2) and (nl_info_comp.n gt (max(nl_info_comp.n)/2 < 15))))
	tmp=fltarr(n_elements(coords)/2) & conts=intarr(n_elements(cont_info))
	cedpt=replicate({pt:lonarr(max(cont_info.n*3)), sep:0, start:0, fin:0}, n_elements(cont_info))

	for i=0, n_elements(cont_info)-1 do begin
		if (cont_info(i).n/3) lt sep then cedpt(i).sep=cont_info(i).n/3 else cedpt(i).sep=sep
		cedpt(i).start=cont_info(i).offset
		cedpt(i).fin=cont_info(i).offset+cont_info(i).n-1
	endfor

		for i=0, n_elements(nl_info_comp)-1 do begin
	for j=nl_info_comp(i).offset, nl_info_comp(i).offset+nl_info_comp(i).n-1 do begin
		x=nl_xy_comp(0,j) & y=nl_xy_comp(1,j)
		tmp=((coords(0,*)-x)^2+(coords(1,*)-y)^2)^.5
		min=min(tmp,min_sub)
		k=0 & while ((min_sub lt cont_info(k).offset) or (min_sub gt (cont_info(k).offset+cont_info(k).n-1))) and k lt n_elements(cont_info)-1 do k=k+1


		if (k lt (n_elements(cont_info)-1)) or ((k eq (n_elements(cont_info)-1)) and (min_sub le (cont_info(k).offset+cont_info(k).n-1)) and (min_sub ge (cont_info(k).offset))) then begin
		conts(k)=conts(k)+1
		cedpt(k).pt(max(where(cedpt(k).pt))+1)=min_sub
		endif
	endfor
		endfor

	for i=0, n_elements(cedpt)-1 do begin
	cedpt(i).pt=cedpt(i).pt(sort(cedpt(i).pt))
	tmp=where(cedpt(i).pt)
	if tmp(0) ne -1 then begin & cedpt(i).pt(0:n_elements(tmp)-1)=cedpt(i).pt(tmp) & cedpt(i).pt(n_elements(tmp):n_elements(cedpt(i).pt)-1)=0 & endif
	endfor

	nl_info_edit=replicate({offset:0L,n:0L,typ:0},(n_elements(nl_info)*20))
	nl_xy_edit=fltarr(2,n_elements(coords)/2)

	for i=0, n_elements(cedpt)-1 do begin
	tmp=cedpt(i) & edpt, tmp & cedpt(i)=tmp
	endfor

	cropcont, cedpt, cont_info, nl_info_edit
	qtst=where(nl_info_edit.n)
	if qtst(0) ne -1 then nl_info_edit=nl_info_edit(qtst) else begin
	no_line=1
	goto, no_lines
	endelse
	trimcoords, coords, nl_info_edit, nl_xy_edit

	endif

; Crops contour at the found endpoints
	cropcont, edpt, cont_info, nl_info

; Create nl_xy structure to hold just the coordinates used by the contours
	nl_xy=fltarr(2, n_elements(coords)/2)

if loud then begin
;print, nl_info(0)
;wt='' & read, wt
;print, coords(*,nl_info(0).offset:nl_info(0).offset+nl_info(0).n-1)
;wt='' & read, wt
endif

; Trim the coords structure down into nl_xy, which will hold only the coordinates used by nl
	trimcoords, coords, nl_info, nl_xy

; Sets the distance across which to cut isthmuses
	if n_elements(isth_cut) eq 0 then isth_cut=1.5

; Checks for isthmuses, cuts out extraneous points and connects contours
; across gap
	isth, nl_info, nl_xy, nl1, isth_cut=isth_cut, thresh=(thresh*1.5), pbt_thresh=(pbt_thresh*1.5)

if loud then begin
;print, nl_info(0)
;wt='' & read, wt
;print, nl_xy(*,nl_info(0).offset:nl_info(0).offset+nl_info(0).n-1)
;wt='' & read, wt
endif

no_lines:
endif

end