;+ 
; NAME: 
; SWARM_PCP_CALC_NGRAD
; 
; PURPOSE: 
; This routine calculates the spatial gradient of the plasma density, in units of m^(-3)/m. It both calculates a running gradient as well as the one for each edge of the polar cap patches. It calculates the gradient as the slope of a linear least-square fit to the (unfiltered) plasma density data.
; 
; CALLING SEQUENCE: 
; swarm_pcp_calc_ngrad, juls, ndens, v_sc, glats, glons, alts, ps1idx, ps2idx, pe1idx, pe2idx
;
; INPUTS:
; juls: The Julian dates (as provided by IDL's JULDAY function) of the data points from the EFI 
; data files. For n data point inside the file this is a vector of length n.
; ndens: The plasma density provided by the Langmuir probe of the EFI, from the EFI data files. For n data point inside the file this is a vector of length n.
; v_sc: The spacecraft velocity, from the EFI data files. For n data point inside the file this is a vector of diminsion [n,3].
; glats: The geographic latitudes of the satellite, from the EFI data files. For n data point inside the file this is a vector of length n.
; glons: The geographic longitudes of the satellite, from the EFI data files. For n data point inside the file this is a vector of length n.
; alts: The altitudes of the satellite, from the EFI data files. For n data point inside the file this is a vector of length n.
; ps1idx: An array of indices which mark the beginning of the patch proper of all patches found.
; ps2idx: An array of indices which mark the beginning of the patch proper of all patches found.
; pe1idx: An array of indices which mark the end of the second edge of all patches found.
; pe2idx: An array of indices which mark the end of the patch proper of all patches found.
;
; OUTPUTS:
; None.
;
; OPTIONAL INPUTS:
; None.
;
; OPTIONAL OUTPUTS:
; None.
;
; KEYWORD PARAMETERS:
; ndens_gradient: Set this to a named variable that will upon completion contain the plasma density gradient calculated as the slope of a linear least-square fit to the (unfiltered) plasma density data. The window width over which the data is fitted is 27 points.
; error_ndens_gradient: Set this to a named variable that will upon completion contain the 1-sigma uncertainty estimate of the slope of a linear least-square fit to the (unfiltered) plasma density data. The error in the measurements (the plasma density) is assumed to be 10%.
; ndens_edge_gradient: Set this to a named variable that will upon completion contain the plasma density gradient calculated as the slope of a linear least-square fit to the (unfiltered) plasma density data. The window width over which the data is fitted is the length of the polar cap patch edge.
; error_ndens_edge_gradient: Set this to a named variable that will upon completion contain the 1-sigma uncertainty estimate of the slope of a linear least-square fit to the (unfiltered) plasma density data on the polar cap patch edges. The error in the measurements (the plasma density) is assumed to be 10%.
;	PGRADLEN: Set this to a named variable that will upon completion contain the number of points over which the gradient is determined -  this will always be returned as 27.
; 
; EXAMPLE: 
; ; read the electric field data, i.e., the data from both the Thermal Ion Imager and the 
; ; Langmuir Probe
; swarm_efi_read, adate, probe,	$
; 	juls=juls, glats=glats, glons=glons, alts=alts, $
; 	ndens=ndens, v_sc=v_sc, v_ion=v_ion, te=te, $
; 	preliminary=preliminary
; 
; ; Calculate the density gradient
; swarm_pcp_calc_ngrad, juls, ndens, v_sc, glats, glons, alts, $
; 	ps1idx, ps2idx, pe1idx, pe2idx, $
; 	ndens_gradient=ndens_gradient, error_ndens_gradient=error_ndens_gradient, $
; 	ndens_edge_gradient=ndens_edge_gradient, error_ndens_edge_gradient=error_ndens_edge_gradient, $
; 	PGRADLEN=PGRADLEN
; 
; COPYRIGHT:
; This file was developed as part of the Swarm + Innovation: Polar Cap Products project. 
; Please visit http://www.mn.uio.no/fysikk/english/research/projects/swarm/ for more information.
; PCP (Polar Cap Products) is a project funded by the European Space Agency 
; (Contract No. 4000114121/15/NL/MP) in the framework of the STSE (Support To Science Element) 
; Swarm + Innovation Program.
; If you have any questions, please contact Lasse Clausen (lasse.clausen@fys.uio.no) or Wojciech 
; Miloch (w.j.miloch@fys.uio.no).
;-
pro swarm_pcp_calc_ngrad, juls, ndens, v_sc, glats, glons, alts, $
	ps1idx, ps2idx, pe1idx, pe2idx, $
	ndens_gradient=ndens_gradient, error_ndens_gradient=error_ndens_gradient, $
	ndens_edge_gradient=ndens_edge_gradient, error_ndens_edge_gradient=error_ndens_edge_gradient, $
	PGRADLEN=PGRADLEN

	PGRADLEN = 27
	pgrad2 = ( PGRADLEN - 1 )/2

	if v_sc[0] eq 0 then begin
		pos = pol2ecef( glats, glons, !re+alts )
		vx = deriv( juls, reform(pos[*,0]) )/86400.d
		vy = deriv( juls, reform(pos[*,1]) )/86400.d
		vz = deriv( juls, reform(pos[*,2]) )/86400.d
		v_sc = [ [vx], [vy], [vz] ]*1e3
	endif

	dsecs = deriv( juls )*86400.d
	vt = sqrt( total( v_sc^2, 2 ) )

	np = n_elements(ndens)
	ndens_gradient = fltarr( np )
	error_ndens_gradient = fltarr( np )
	ndens_edge_gradient = fltarr( np )
	error_ndens_edge_gradient = fltarr( np )
	for i=pgrad2, np-pgrad2-1L do begin
		xx = juls[i-pgrad2:i+pgrad2]
		yy = ndens[i-pgrad2:i+pgrad2]*1e6
		tmp = linfit( xx, yy, measure_error=0.1*yy, sigma=sigma )
		ndens_gradient[i] = tmp[1]/86400.d/mean( vt[i-pgrad2:i+pgrad2] )
		error_ndens_gradient[i] = sigma[1]/86400.d/mean( vt[i-pgrad2:i+pgrad2] )
	endfor

	n_p = n_elements(ps1idx)
	for p=0, n_p-1 do begin

;		if p eq 0 then si = 0 else si = pe1idx[p-1]+1
;		if p eq n_p-1 then fi = np-1 else fi = ps1idx[p]-1
;		if fi-si ge 0 then begin
;			ndens_gradient[si:fi] = 0.
;			error_ndens_gradient[si:fi] = 0.
;		endif

		xx = juls[ps1idx[p]:ps2idx[p]]
		yy = ndens[ps1idx[p]:ps2idx[p]]*1e6
		tmp = linfit( xx, yy, measure_error=0.1*yy, sigma=sigma )
		ndens_edge_gradient[ ps1idx[p]:ps2idx[p] ] = tmp[1]/86400.d/mean( vt[ps1idx[p]:ps2idx[p]] )
		error_ndens_edge_gradient[ ps1idx[p]:ps2idx[p] ] = sigma[1]/86400.d/mean( vt[ps1idx[p]:ps2idx[p]] )

		xx = juls[pe2idx[p]:pe1idx[p]]
		yy = ndens[pe2idx[p]:pe1idx[p]]*1e6
		tmp = linfit( xx, yy, measure_error=0.1*yy, sigma=sigma )
		ndens_edge_gradient[ pe2idx[p]:pe1idx[p] ] = tmp[1]/86400.d/mean( vt[pe2idx[p]:pe1idx[p]] )
		error_ndens_edge_gradient[ pe2idx[p]:pe1idx[p] ] = sigma[1]/86400.d/mean( vt[pe2idx[p]:pe1idx[p]] )

	endfor

end