;Derived functions/procedures for the W05SC model
;Uses analytic SCHA functions
;This is a shorter version that does not include geographic coordinate conversions
;**********************************************
Function dPm_n,m,n,cth,dbl=dbl ;Calculates the derivative of Pm_n
;allow for n to be an array
sth=sqrt(1.D -cth^2)
xn=n*(n+1)
x=(1.D -cth)/2.D
nresult=n_elements(cth)>n_elements(n)
IF m eq 0 Then Begin
  A =.5*sth
  dA=replicate(0.d,nresult)
endif else begin
  IF m eq 1 Then Begin
    A=sth *.5*sth
    dA=cth
  endif Else Begin
    A=sth^m *.5*sth
    dA=m*sth^(m-1) * cth
  Endelse
endelse
if n_elements(cth) LT n_elements(n) then Begin
  A=replicate(A[0],nresult)
  dA=replicate(dA[0],nresult)
  x=replicate(x[0],nresult)
endif
result=dA
k=1
REPEAT Begin
  f=((k+m-1)*(k+m)-xn) / (k*(k+m))
  fx=f*x
  IF k EQ 1 Then A *=f Else A *= fx
  dA *= fx
  delta=dA+A*k
  result += delta
  k++
EndREP UNTIL MAX(ABS(delta)/(ABS(result)>1.e-6)) LT 1.E-6
if keyword_set(dbl) then return,KM_n(m,n)*result else Return,float(KM_n(m,n)*result)
end
;**********************************************
Function SCdPlm,i,colat,nlm
;return derivative of Spherical Cap Harmonic Associated Legendre values,
;given colat values and index i into array of L and M values
;reduce computational overhead by computing all spherical cap harmonic Plm values
;in this routine, so that if it is called again with identical input, then they
;do not need to be recomputed again.
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize
common saveSCplm,prevth0,ColatTable,PlmTable,nlms
common saveSCdplm,myth0,dPlmTable
COMMON SetBoundaryData,BndyFitR
th0=BndyFitR[0]
if n_elements(prevth0) EQ 0 || th0 NE prevth0 Then junk=SCplm(i,0.) ;force creation of new ColatTable
if n_elements(myth0) EQ 0 || th0 NE myth0 then begin
  TableSize=n_elements(ColatTable)
  dPlmTable=fltarr(TableSize,csize)
  myth0=th0
  cth=COS(ColatTable*(!DPI/180.))
  For j=0,csize-1 Do Begin
    m=Ms[j]
    dPlmTable[0,j]=dPm_n(m,nlms[j],cth)
    if m ne 0 && AB[j] then begin ;AB=1 always preceeds 0 with same l & m values
      dPlmTable[0,j+1]=dPlmTable[*,j]
      j++
    endif
  endfor
endif
nlm=nlms[i]
return,interpol(dPlmTable[*,i],ColatTable,colat,/QUAD)
end
;**********************************************
FUNCTION MPFAC,lat,mlt,R,ALT=alt,INSIDE=inside,OUTSIDE=outside,FILL_VALUE=fill
;Returns the FAC from the Magnetic Potential model
;Evaluate the Laplacian (Grad^2) of the potential, at radius R
;Uses a spherical harmonic recursion formula
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc
COMMON SetBoundaryData,BndyFitR

IF n_elements(alt) NE 0 then R=6371.2+alt
IF n_elements(R) EQ 0 Then R=6371.2+110. ;km radius

precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phir
if arg_present(outside) then outside=inside EQ 0
if n_elements(fill) eq 0 then fill=0.
if ngood EQ 0 Then return,make_array(nresult,VALUE=fill)

  z=fltarr(ngood)

  phim=phir # replicate(1,maxm) * ((indgen(maxm)+1) ## replicate(1,n_elements(phir)))
  cospm=cos(phim)
  sinpm=sin(phim)

  For j=0,csize-1 Do Begin
    if Ls[j] ge 11 then break ;only go up to L=10 for FAC
    m=Ms[j]
    if AB[j] then begin ;AB=1 always preceeds 0 with same l & m values
      Plm=SCPlm(j,colats,nlm)
      Plm *= nlm*(nlm+1) ;recursion formula for the 2d Laplacian operation
      if m eq 0 then Z-= Plm*Bsphc[j] $
      else begin
        Z-= Plm*(Bsphc[j]*cospm[*,m-1]+Bsphc[j+1]*sinpm[*,m-1])
        j++
      endelse
    endif
  endfor

cfactor= -1.e5/(4.*!PI*R^2)   ;convert to uA/m2
;sign has convention that negative is upward
Z *= cfactor

if keyword_set(Geographic) then Z *=ascale

if ngood EQ nresult Then RETURN,Z
;otherwise:
result=make_array(nresult,VALUE=fill)
result[goodpts]=Z
 
;print,format="('mpfac: lat=',f8.3,' mlt=',f8.3,' fac=',e12.4)",$
;  lat,mlt,result

Return,result
END
;*************************************************************************
FUNCTION EpotFAC,lat,mlt,R,Pedersen=pedersen,INSIDE=inside,OUTSIDE=outside,$
  FILL_VALUE=fill
;Evaluate the Laplacian (Grad^2) of the potential, at radius R
;and returns the FAC that would exist with a uniform Pedersen conductivity
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc
COMMON SetBoundaryData,BndyFitR

IF n_elements(pedersen) EQ 0 Then pedersen=7. ;mhos
IF n_elements(R) EQ 0 Then R=6371.2+110. ;km radius

precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phir
if arg_present(outside) then outside=inside EQ 0
if n_elements(fill) eq 0 then fill=0.
if ngood EQ 0 Then return,make_array(nresult,VALUE=fill)

z=fltarr(ngood)

  phim=phir # replicate(1,maxm) * ((indgen(maxm)+1) ## replicate(1,n_elements(phir)))
  cospm=cos(phim)
  sinpm=sin(phim)

  For j=0,csize-1 Do Begin
    if Ls[j] ge 11 then break ;only go up to L=10 for FAC
    m=Ms[j]
    if AB[j] then begin ;AB=1 always preceeds 0 with same l & m values
      Plm=SCPlm(j,colats,nlm)
      Plm *= nlm*(nlm+1)
      if m eq 0 then Z-= Plm*Esphc[j] $
      else begin
        Z-= Plm*(Esphc[j]*cospm[*,m-1]+Esphc[j+1]*sinpm[*,m-1])
        j++
      endelse
    endif
  endfor


cfactor= -1000.*pedersen/(R^2)   ;convert to uA/m2
;sign has convention that negative is upward
if ngood EQ nresult Then RETURN,Z*cfactor
;otherwise:
result=make_array(nresult,VALUE=fill)
result[goodpts]=Z*cfactor
Return,result
END
;**********************************************
PRO EMFIELDS,lat,mlt,ALT=alt,SOUTHPOLE=southpole,$
  BNORTH=bnorth,BEAST=Beast,ENORTH=enorth,EEAST=eeast,JNORTH=jnorth,JEAST=jeast,$
  INSIDE=inside,OUTSIDE=outside,FILL_VALUE=fill
;Computes Electric field, Magnetic field, and/or the ionospheric currents,
;  from the electric and magnetic potentials,
;  at given latitude/MLT and specified altitude.
;BNORTH & BEAST = components of the magnetic perturbation due to FAC, above
;  the ionosphere, as would be measured with a satellite.
;ENORTH & EEAST = components of the electric field at & above ionosphere.
;JNORTH & JEAST = components of the total Pedersen current in the ionosphere.
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc

Re=6371.2
default,alt,110.
R=Re+alt
default,southpole,0

precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phir,ANGLE=angle
if arg_present(outside) then outside=inside EQ 0
if n_elements(fill) eq 0 then fill=0.
if ngood EQ 0 Then Begin ;all lat-mlt points are outside of the model's boundary circle
  zeros=make_array(nresult,VALUE=fill)
  if arg_present(Bnorth)then Bnorth=zeros
  if arg_present(Jnorth)then Jnorth=zeros
  if arg_present(Enorth)then Enorth=zeros
  if arg_present(Beast) then Beast=zeros
  if arg_present(Jeast) then Jeast=zeros
  if arg_present(Eeast) then Eeast=zeros
  return
endif

if ngood NE nresult then zeros=make_array(nresult,VALUE=fill)
;if lat and mlt are not sized the same, then get dimensions of the larger array,
;  which assumes that the other one is a scalar.
if n_elements(lat) GE n_elements(mlt) Then dim=SIZE(lat,/DIMENSION) Else dim=SIZE(mlt,/DIMENSION)

Sina=SIN(angle)
Cosa=COS(angle)

SineColats=SIN(colats*!DTOR) >1.e-8

phim=phir # replicate(1,maxm) * ((indgen(maxm)+1) ## replicate(1,n_elements(phir)))
cospm=cos(phim)
sinpm=sin(phim)

if arg_present(Bnorth) || arg_present(Beast) then begin
  BN=fltarr(ngood)
  BE=fltarr(ngood)
  domagnetic=1
endif else domagnetic=0
IF arg_present(Jnorth)|| arg_present(Jeast) then begin
  JE=fltarr(ngood)
  JN=fltarr(ngood)
  doJperp=1
endif else doJperp=0
if arg_present(Enorth)|| arg_present(Eeast) then begin
  EE=fltarr(ngood)
  EN=fltarr(ngood)
  doElectric=1
endif else doElectric=0


For j=0,csize-1 Do Begin
  l=Ls[j] & m=Ms[j]
  IF l NE 0 Then Begin
    IF m EQ 0 Then Begin
      dPlms=SCdPlm(j,colats)
      If domagnetic then BE += Bsphc[j]*dPlms
      if doJperp then JN += Bsphc[j]*dPlms
      if doelectric then EN += Esphc[j]*dPlms
    endif else begin
      dPlms=SCdPlm(j,colats)
      Plms=SCPlm(j,colats)
      mm=m-1
      if domagnetic then begin
        A1=Bsphc[j]
        A2=Bsphc[j+1]
        BE += (A1*cospm[*,mm] + A2*sinpm[*,mm])*dPlms
        BN += (A2*cospm[*,mm] - A1*sinpm[*,mm])*m*Plms
      endif
      if doJperp then begin
        A1=Bsphc[j]
        A2=Bsphc[j+1]
        JN += (A1*cospm[*,mm] + A2*sinpm[*,mm])*dPlms
        JE += (A2*cospm[*,mm] - A1*sinpm[*,mm])*m*Plms
      endif
      if doelectric then begin
        A1=Esphc[j]
        A2=Esphc[j+1]
        EN += (A1*cospm[*,mm] + A2*sinpm[*,mm])*dPlms
        EE += (A2*cospm[*,mm] - A1*sinpm[*,mm])*m*Plms
      endif
      j ++
    endelse
  Endif
Endfor

if domagnetic then begin
  ;get North and East magnetic field
  f=10000./R
  BN *= f/SineColats
  IF southpole THEN f *= -1.
  BE *= f

  ;now the North and East components need to be rotated from the offset pole system
  ;back into the original CGM coordinate system.
  if ngood EQ nresult Then Begin
    Bnorth=BN*Cosa+BE*Sina
    Beast=BE*Cosa-BN*Sina
  endif else begin
    Bnorth=zeros
    Beast=zeros
    Bnorth[goodpts]=BN*Cosa+BE*Sina
    Beast[goodpts]=BE*Cosa-BN*Sina
  endelse

  if n_elements(dim) GT 1 Then Begin
    Bnorth=reform(Bnorth,dim,/overwrite)
    Beast=reform(Beast,dim,/overwrite)
  endif
Endif ;domagnetic

if doJperp then begin
  ;get North and East current
  f= -1.e5/(4.*!PI*R)   ;converts cTm/km to A/km
  JE *= f/SineColats
  IF southpole EQ 0 THEN f= -1.*f
  JN= f*temporary(JN)

  ;now the North and East components need to be rotated back into the original
  ;CGM coordinate system.
  if ngood EQ nresult Then Begin
    Jnorth=JN*Cosa+JE*Sina
    Jeast=JE*Cosa-JN*Sina
  endif else begin
    Jnorth=zeros
    Jeast=zeros
    Jnorth[goodpts]=JN*Cosa+JE*Sina
    Jeast[goodpts]=JE*Cosa-JN*Sina
  endelse

  if n_elements(dim) GT 1 Then Begin
    Jeast=reform(Jeast,dim,/overwrite)
    Jnorth=reform(Jnorth,dim,/overwrite)
  endif
endif ;doJperp

if doElectric then begin
  ;get North and East electric field
  f= -1000./R   ;converts to mV/m
  EE *= f/SineColats
  IF southpole EQ 0 THEN f= -1.*f
  EN *= f

  ;now the North and East components need to be rotated back into the original
  ;CGM coordinate system.
  if ngood EQ nresult Then Begin
    Enorth=EN*Cosa+EE*Sina
    Eeast=EE*Cosa-EN*Sina
  endif else begin
    Enorth=zeros
    Eeast=temporary(zeros)
    Enorth[goodpts]=EN*Cosa+EE*Sina
    Eeast[goodpts]=EE*Cosa-EN*Sina
  endelse

  if n_elements(dim) GT 1 Then Begin
    Eeast=reform(Eeast,dim,/overwrite)
    Enorth=reform(Enorth,dim,/overwrite)
  endif
endif ;doElectric
Return
END
;**********************************************
PRO BFIELD,lat,mlt,Bnorth,Beast,MAGNITUDE=magnitude,_REF_EXTRA=extra
EMFIELDS,lat,mlt,BNORTH=bnorth,BEAST=beast,_EXTRA=extra
if keyword_set(magnitude)then Bnorth=sqrt(Beast^2+temporary(Bnorth)^2)
RETURN
END
;**********************************************
PRO Jperp,lat,mlt,Jnorth,Jeast,MAGNITUDE=magnitude,_REF_EXTRA=extra
EMFIELDS,lat,mlt,JNORTH=jnorth,JEAST=jeast,_EXTRA=extra
if keyword_set(magnitude)then Jnorth=sqrt(Jeast^2+temporary(Jnorth)^2)
Return
END
;**********************************************************************
PRO Efield,lat,mlt,Enorth,Eeast,MAGNITUDE=magnitude,_REF_EXTRA=extra
EMFIELDS,lat,mlt,ENORTH=enorth,EEAST=eeast,_EXTRA=extra
if keyword_set(magnitude)then Enorth=sqrt(Eeast^2+temporary(Enorth)^2)
Return
END
;**********************************************
function JouleHeat,Lat,Mlt,FILL=fill,OUTSIDE=outside, _REF_EXTRA=extra
;Returns the Joule Heating in mW/m2

EMFIELDS,lat,mlt,ENORTH=enorth,EEAST=eeast,Jnorth=Jnorth,Jeast=Jeast,$
    OUTSIDE=outside,FILL=0.,_EXTRA=extra

result=1.e-3*(Jnorth*Enorth + Jeast*Eeast) > 0.
;almost always zero, the plots are smoother if the few points with
;low-magnitude negative values are removed.


if n_elements(fill) eq 0 then fill=0.
if fill eq 0. then return,result
;otherwise, set outside points to provided fill value
fillpoints=where(outside,nfill)
if nfill eq 0 then return,result
result[fillpoints]=fill
return,result
End
;**********************************************
Pro W05derived,version
version='V4bsc'
end