; 2005 Version of the electric and magnetic potential (FAC) models
; by Dan Weimer.  Uses Spherical Cap Harmonic Analysis (SCHA) functions.
; Model description is in:
; Weimer, D. R., Predicting Surface Geomagnetic Variations Using Ionospheric
;    Electrodynamic Models, J. Geophys. Res., 110, A12307, doi:10.1029/
;    2005JA011270, 2005.
; Some information about the model (such as outer boundary calculation)
; is also in the earlier paper:
; Weimer, D. R. (2005), Improved ionospheric electrodynamic models and
;    application to calculating Joule heating rates,  J. Geophys. Res., 110,
;    A05306, doi:10.1029/2004JA010884.
;
; For information about the SCHA, see the paper:
;Haines, G. V., Spherical cap harmonic analysis, J. Geophys. Res., 90, B3,
;  2583, 1985.  (Note that this is in JGR-B, "Solid Earth", rather than JGR-A)
;
;All data files are in the IDL SAVE/RESTORE XDR format.  For conversion to
;FORTRAN it will be neccessary to use IDL to RESTORE each file, then write out
;all varriables using WRITEU, in a FORTRAN and machine compatible format,
;perhaps with big- and little-endian versions.
;**************************************************************************
;@DEFAULT.pro
@PATH_SEP.pro ;in IDL Library
;@SOURCEPATH.pro ;function that tells IDL what directory contains the
; .pro source code or .sav executable that is being run.  Allows for the data
; files to be located with the source code, while having a different default
; directory.
@UNIQ.pro ;in IDL Library
;*****************************************************************************
; 3/26/08 btf: my interpretation of pro default (did not receive DEFAULT.pro 
;              from Weimer, which was included above)
pro default,var,value
  if (n_elements(var) eq 0) then var = value
end
;*****************************************************************************
; 3/26/08 btf: my interpretation of func sourcepath (did not receive 
;              SOURCEPATH.pro from Weimer, which was included above)
; restore,sourcepath('nkmlookup',/functions)+'SCHAtable.xdr'
;
function sourcepath,type,functions=functions
  return,'./' ; assume all source is in cwd
end
;*****************************************************************************
Function Km_n,m,n
;A normalization function used by the SCHA routines.  See Haines.
;If n is an array, then returns an array.
IF m eq 0 Then Return,1.D
good=WHERE(n GE m,ngood)
if ngood eq n_elements(n) Then $
return,SQRT( 2.D * exp(lngamma(double(n+m+1))-lngamma(double(n-m+1)))) / (2.D ^m * factorial(m) )
;or return,SQRT( 2.D *factorial(n+m)/factorial(n-m) ) / (2.D ^m * factorial(m) )
sz=size(n)
if sz[0] eq 0 then result=0.d Else result=make_array(SIZE=sz,/double)
result[good]=$
  SQRT( 2.D * exp(lngamma(double(n[good]+m+1))-lngamma(double(n[good]-m+1)))) / (2.D ^m * factorial(m) )
return,result
End
;*****************************************************************************
Function Pm_n,m,n,cth,dbl=dbl
;Another SCHA function, returns the SCHA version of the associated Legendre
;Polynomial, Pmn
;If n is an array, then returns an array.
nresult=n_elements(cth)>n_elements(n)
IF m EQ 0 Then A=replicate(1.d,nresult) Else A=sqrt(1.D -cth^2)^m
xn=n*(n+1)
x=(1.D -cth)/2.D
if n_elements(cth) LT n_elements(n) then Begin
  A=replicate(A[0],nresult)
  x=replicate(x[0],nresult)
endif

;print,format="('Pm_n: a=',/,(6e12.4))",a
;print,format="('Pm_n: xn=',e12.4)",xn
;print,format="('Pm_n: x=',/,(6e12.4))",x

result=A
k=1
REPEAT Begin
  A *= x*((k+m-1)*(k+m)-xn)/(k*(k+m))

; print,format="('Pm_n: k=',i3,' a=',/,(6e12.4))",k,a

  result += A

; print,format="('Pm_n: k=',i3,' result=',/,(6e12.4))",k,result

  k++

; print,format="('Pm_n: k=',i3,' abs(a)=',/,(6e12.4))",k,abs(a)
; print,format="('Pm_n: k=',i3,' abs(table)=',/,(6e12.4))",k,abs(result)
  
  tmp = abs(a)/(abs(result)>1.e-6)
; print,format="('Pm_n: k=',i3,' tmp=',/,(6e12.4))",k,tmp
; print,format="('Pm_n: k=',i3,' max(tmp)=',e12.4)",k,max(tmp)

EndREP UNTIL MAX(ABS(A)/(ABS(result)>1.e-6)) LT 1.E-6
kmn = Km_n(m,n)

;print,format="('Pm_n: kmn=',e12.4,' result=',/,(6e12.4))",kmn,result

if keyword_set(dbl) then return,KM_n(m,n)*result else Return,float(KM_n(m,n)*result)
end
;*****************************************************************************
Function nkmlookup,k,m,th0
;Given the size of a spherical cap, defined by the polar cap angle, th0, and
;also the values of integers k and m, returns the value of n, a real number (see Haines).
;It uses interpolation from a lookup table that had been precomputed, in order
;to reduce the computation time.
common nkmTabDat,th0s,allnkm,maxk,maxm
if th0 EQ 90. then return,k
;uses file created with "Allnksearch" procedure.
;print,'nkmlooup: th0s=',th0s
if n_elements(th0s) eq 0 Then begin
  path = sourcepath('nkmlookup',/functions)+'SCHAtable.xdr'
  restore,sourcepath('nkmlookup',/functions)+'SCHAtable.xdr'
endif
if k gt maxk then begin
  print,'Warning: nkmlookup called with k out of range of table'
  result=interpol(allnkm[maxk,m,*],th0s,th0,/quadratic)
  return,(result[0] > k)
endif
if m gt maxm then begin
  print,'Warning: nkmlookup called with m out of range of table'
  result=interpol(allnkm[k,maxm,*],th0s,th0,/quadratic)
  return,(result[0] > k)
endif
test=where(th0 LT th0s[0],ntest)
if ntest GT 0 then print,'Warning: nkmlookup called with th0 less than lowest in table.'

;print,format=$
;  "('nkmlookup call interpol: k=',i3,' m=',i3,' th0=',e12.4,' allnkm=',/,(6e12.4))",$
;  k,m,th0,allnkm[k,m,*]

result=interpol(allnkm[k,m,*],th0s,th0,/quadratic)
if (size(th0))[0] eq 0 then result=result[0]

return_val = (result > k)
;help,return_val
;print,format="('nkmlookup: k=',i3,' m=',i3,' return_val=',(f10.4))",$
;  k,m,return_val

return,(result > k)
end
;*********************************************************************************
PRO SetBoundary,angle,Bt,Tilt,SWVel,SWDen,ALindex,YZ=yz,PATH=mypath
;Sets the coefficients that define the low-latitude boundary model,
;  given the IMF and solar wind values.
;if keyword YZ is set, then the first two parameters are By and Bz,
; in place of the 'angle' and BT parameters.

COMMON SaveW05bBndyCoef,BndyA,BndyB,ex
COMMON SetBoundaryData,BndyFitR,Tmat,TTmat
common geoboundarydata,needupdate,bndyglats,bndyglons ;used by function geoboundarylat
needupdate=1
if n_elements(BndyB) eq 0 then Begin
  ;has BndyA,BndyB
  if ~keyword_set(mypath) then mypath=sourcepath('SetBoundary')
  file=mypath+'W05scBndy.xdr'
  RESTORE,file
  ;calculate the transformation matrix to the coordinate system of the offset pole
  xc=4.2D0
  theta=xc*(!DPI/180.)
  CT=cos(theta)
  ST=sin(theta)
  TMat=$ ;Transformation matrix
  [[ CT ,0.d , ST ],$
   [0.d ,1.d ,0.d ],$
   [-ST ,0.d , CT ]]
  TTMat=TRANSPOSE(TMat)
endif

if keyword_set(yz) then begin
  default,angle,0.
  default,Bt,-5.
  p1=angle
  p2=bt
  Bt=sqrt(p1^2+p2^2)
  angle=ATAN(p1,p2)*!RADEG
  IF angle LT 0. THEN angle += 360.
endif else begin
  default,angle,180.
  default,Bt,5.
endelse
default,Tilt,0.
default,SWVel,450.
default,SWDen,8.
swp=SWDen*SWVel^2 *1.6726e-6  ;pressure
tilt2=tilt^2

usesAL= n_elements(ALindex) GT 0
cosa=cos(angle*!DTOR)

btx= 1.-exp(-bt*ex[0])
if bt gt 1 then btx *= bt^ex[1] else $
;if bt lt 1. then $
  cosa=1.+bt*(cosa-1.) ;remove angle dependency for IMF under 1 nT

if usesAL then Begin
  x=[ 1., cosa, btx, btx*cosa, SWVel, swp , ALindex ]
  C=BndyB
endif Else Begin
  x=[ 1., cosa, btx, btx*cosa, SWVel, swp ]
  C=BndyA
endelse
BndyFitR=TOTAL(x*C)

if keyword_set(yz) then begin
 ;return inputs to original values
  angle=p1
  Bt=p2
endif
RETURN
END
;**********************************************************************************
PRO DoRotation,latin,lonin,latout,lonout,nresult,REVERSED=reversed
;uses transformation matrices saved in common block, to convert between
;the given geomagnetic latatud/longitude, and the coordinate system that is
;used within the model,that is offset from the pole.
COMMON SetBoundaryData,BndyFitR,Tmat,TTmat
;Rotate Lat/Lon spherical coordinates with the transformation given by saved matrix
;The coordinates are assumed to be on a sphere of Radius=1.
;Uses cartesian coordinates as an intermediate step.

  DDTOR=!DPI/180.
  LATR=latin*DDTOR
  LONR=lonin*DDTOR
  STC = SIN(LATR)
  CTC = COS(LATR)
  SF = SIN(LONR)
  CF = COS(LONR)
  A=CTC*CF ;if either latin or lonin was scalar and other was vector
  B=CTC*SF ;then the result of the multiplication has same size as the vector
  nresult=N_ELEMENTS(A)
  ;Make sure that all are same size and reformed to one-dimensional
  A=reform(A,nresult,/overwrite)
  B=reform(B,nresult,/overwrite)
  if n_elements(STC) lt nresult then STC=REPLICATE(STC[0],nresult) $
  else STC=reform(stc,nresult,/overwrite)
  if keyword_set(reversed) then TM=Tmat Else TM=TTmat
  ;;The reform to 1D array is required for this:
  ;Pos= TM ## [[CTC*CF],[CTC*SF],[STC]]
  Pos= TM ## [[A],[B],[STC]]
  Latout=reform(FLOAT(ASIN(Pos[*,2])/DDTOR))
  Lonout=reform(FLOAT(ATAN(Pos[*,1],Pos[*,0])/DDTOR))
; print,format="('DoRotation: latin,lonin=',2f9.4,' latout,lonout=',2f9.4)",$
;   latin,lonin,latout,lonout
RETURN
END
;*********************************************************************************
FUNCTION BoundaryLat,mlt
;Given the MLT, returns the (corrected) geomagnetic latitude of the model's
;boundary, in degrees.
COMMON SetBoundaryData,BndyFitR
Rcirc=BndyFitR[0] ;the radius of the boundary in degrees
x0=4.2 ;& y0=0.
phi=mlt*!PI/12.
b=x0*cos(phi);+y0*sin(phi)
;r= b + SQRT(b^2 + Rcirc^2 -x0^2 -y0^2 )
r= b + SQRT(b^2 + Rcirc^2 -x0^2 ) ;y0=0 ;the co-latitude
RETURN,90.-r
END
;*********************************************************************************
PRO READMODEL,junk
;not used now, kept in for compatibility with older versions of the model
RETURN
END
;*********************************************************************************
PRO SetModel,angle,Bt,Tilt,SWVel,SWDen,ALindex,RESET=reset,VERSION=version,YZ=yz,PATH=mypath
;
; Calculate the complete set of the models' SCHA coeficients,
; given an aribitrary IMF angle (degrees from northward toward +Y),
; magnitude Bt (nT), dipole tilt angle (degrees), and
; solar wind velocity (km/sec).
;
; If keyword YZ is set, then the first two parameters are By and Bz,
; in place of the 'angle' and BT parameters.
;
; This procedure creates the array Esphc, in COMMON SetW05sc
; The reset option was only used during development, with a comparison of different
;  model versions that shared the same code.
COMMON SaveW05scEpot,EX,d1,d2,SCHFits,ALSCHFits
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc

if n_elements(SCHFits) eq 0 or keyword_set(Reset) then Begin
  if ~keyword_set(mypath) then mypath=sourcepath('SetModel')
  file=mypath+'W05SCEpot.xdr'
  RESTORE,file
endif
;This file contains:
;CSIZE=28  ;defines size of arrays
;MAXL=12   ;largest L value used in the Plm polynomials
;MAXM=2    ;larmest M value used in the Plm polynomials
;D1=15 ;number of elements in normal model
;D2=18 ;number of elements in optional model using AL index
;EX=Fltarr(2) ;constants for an EXponential equation for saturation effect
;AB=Bytarr(CSIZE) ;denotes if coef. is for a sine or cosine term
;LS=Intarr(CSIZE) ;array of L values
;MS=Intarr(CSIZE) ;array of M values
;Only terms with Ls+Ms=odd number are used in this series
;SCHFITS=Dblarr(d1, CSIZE)
;ALSCHFITS=Dblarr(d2, CSIZE)

version='W05SC'

if keyword_set(yz) then begin
  default,angle,0.
  default,Bt,-5.
  p1=angle
  p2=bt
  Bt=sqrt(p1^2+p2^2)
  angle=ATAN(p1,p2)*!RADEG
  IF angle LT 0. THEN angle += 360.
endif else begin
  default,angle,180.
  default,Bt,5.
endelse
default,Tilt,0.
default,SWVel,450.
default,SWDen,9.

SetBoundary,angle,Bt,Tilt,SWVel,SWDen,ALindex,PATH=mypath

stilt=sin(tilt*!dtor)
stilt2=stilt^2
sw=Bt*swvel/1000.
swe= (1.-exp(-sw*ex[1]))*sw^ex[0]
C0=1.
SWP=SWvel^2 * SWDen*1.6726E-6

Rang=angle*!DTOR
cosa=cos(Rang)
sina=sin(Rang)
cos2a=cos(2*Rang)
sin2a=sin(2*Rang)
if Bt LT 1. then begin ;remove angle dependency for IMF under 1 nT
  cosa= -1.+bt*(cosa+1.)
  cos2a=1.+bt*(cos2a-1.)
  sina= bt*sina
  sin2a=bt*sin2a
endif

;print,'SetModel: swe=',swe,' stilt=',stilt,' stilt2=',stilt2,' cosa=',cosa,$
;  ' sina=',sina,' cos2a=',cos2a,' sin2a=',sin2a

;The optional use of the AL Index does not need to be included in the FORTRAN version.
if n_elements(ALindex) GT 0 then Begin
  al=ALindex[0] ;should be a scalar
  cfits=ALSCHFits  ;ALSCHFits=FLTARR(d2,csize)
  A=[C0, swe, stilt, stilt2, swp, al , $
     swe*cosa, stilt*cosa, stilt2*cosa, swp*cosa, al*cosa,$
     swe*sina, stilt*sina, stilt2*sina, swp*sina, al*sina,$
     swe*cos2a,$
     swe*sin2a ]
  d=d2
endif else begin
  cfits=SCHFits  ;;SCHFits=FLTARR(d1,csize)
  ;;A=[C0, swe, swe*cosa, swe*sina, swe*cos2a, swe*sin2a ]
  A=[C0, swe, stilt,stilt2,swp,$
     swe*cosa, stilt*cosa, stilt2*cosa, swp*cosa,  $
     swe*sina, stilt*sina, stilt2*sina, swp*sina,  $
     swe*cos2a,$
     swe*sin2a ]
  d=d1
endelse

Esphc=fltarr(csize)
for i=0,d-1 do Esphc += reform(cfits[i,*])*A[i]

;print,format="('SetModel: Esphc=',/,(6e12.4))",Esphc

;Added in final release to set up both models at same time
;Not needed here if not using the FAC model.
SetMPfac,angle,Bt,Tilt,SWVel,SWDen,ALindex,PATH=mypath,/SKIPBNDYSET,RESET=reset
;This sets the array Bsphc

if keyword_set(yz) then begin
 ;return inputs to original values
  angle=p1
  Bt=p2
endif

RETURN
END
;*********************************************************************************
Pro precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phiangle,$
  GOODLATS=goodlats,GOODMLTS=goodmlts,RADII=radii,ANGLE=angle
;This procedure performs some opperations on the input vals. of lat (deg.)
;and MLT (hours), as needed by the electric and magnetic potential functions.
;Some of the results are required by graphic programs, to determine from a given
;set of input values which points are within the models' boundary.
;The outputs are:
;  nresult: the number of elements in the result, since either lat or mlt may be
;    given as scalars or arrays.
;  ngood: the number of valid points, with locations within the models' boundary.
;  goodpts: array of indices to valid points.
;  inside: logical array of size nresult, set to 1 for each point within the boundary.
;  radii:  radii/colatitude within the models'; offset coordinate system.
; The following arrays (or scalar) that are returned have a size = ngood
;  colats: radii[goodpts]
;  phiangle: longitude in RADIANS within the models' offset coordinate system.
;  goodlats: array that contains only the input latitudes that are within the boundary
;  goodmlts: same thing for mlt
;  angle: difference between angles. Used when calculating electric and magnetic field vectors
COMMON SetBoundaryData,BndyFitR
lons=mlt*15.
DoRotation,lat,lons,Tlat,Tlon,nresult
radii=90.-Tlat
inside=radii LE BndyFitR[0] ;Bndycolat
goodpts=where(inside,ngood)
if ngood EQ 0 then return
if arg_present(colats) then colats=radii[goodpts]
if arg_present(phiangle) then phiangle=Tlon[goodpts]*!DTOR
if arg_present(angle) then angle=(Tlon[goodpts]-lons[goodpts])*!DTOR
if arg_present(goodlats) then begin
  if n_elements(lat) lt nresult then goodlats=replicate(lat[0],ngood) $
  else goodlats=lat[goodpts]
endif
if arg_present(goodmlts) then begin
  if n_elements(mlt) lt nresult then goodmlts=replicate(mlt[0],ngood) $
  else goodmlts=mlt[goodpts]
endif
return
end
;*********************************************************************************
Function SCPlm,i,colat,nlm,dbl=dbl
;return 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 SetBoundaryData,BndyFitR

th0=BndyFitR[0]

if n_elements(prevth0) EQ 0 || th0 NE prevth0 Then Begin
  TableSize=3*round(th0)

; print,format="('SCPlm: index=',i3,' colat=',f8.3,' th0=',e12.4,' tablesize=',i3)",$
;   i,colat,th0,tablesize

  PlmTable=dblarr(TableSize,csize)
  nlms=fltarr(csize)
  ColatTable=findgen(TableSize)*(th0/(TableSize-1))
  prevth0=th0
  cth=COS(ColatTable*(!DPI/180.))

; print,format="('SCPlm: colattable=',/,(6f8.3))",colattable
; print,format="('SCPlm: cth=',/,(6e12.4))",cth

  For j=0,csize-1 Do Begin
    l=Ls[j] & m=Ms[j]
    nlms[j]=nkmlookup(l,m,th0)

;   print,format=$
;     "('SCPlm after nkmlookup: j=',i3,' l=',i3,' m=',i3,' nlms[j]=',f8.4)",$
;     j,l,m,nlms[j]

    PlmTable[0,j]=Pm_n(m,nlms[j],cth,/dbl)

;   print,format="('SCPlm: j=',i3,' PlmTable[0,j]=',e12.4)",j,PlmTable[0,j]
;   help,PlmTable[0,j]

    if m ne 0 && AB[j] then begin ;AB=1 always preceeds 0 with same l & m values
      PlmTable[0,j+1]=PlmTable[*,j]
      nlms[j+1]=nlms[j]
      j++
    endif
  endfor
endif
nlm=nlms[i]

;help,PlmTable
;print,format="('SCPlm: PlmTable min,max=',2e12.4)",min(PlmTable),max(PlmTable)

result=interpol(PlmTable[*,i],ColatTable,colat,/QUAD)

;print,format="('SCPlm: i=',i4,' result=',e12.4,' PlmTable[*,i]=',/,(6e12.4))",$
;  i,result,PlmTable[*,i]

;print,format="('SCPlm returning: i=',i4,' result=',e12.4)",i,result

if keyword_set(dbl) then return,result else return,float(result)
end
;*********************************************************************************
FUNCTION EpotVal,lat,mlt,INSIDE=inside,OUTSIDE=outside,FILL_VALUE=fill,BPOT=bpot
;Return the Potential (in kV) at given combination of def. latitude (lat) and MLT,
; in geomagnetic apex coordinates (practically identical to AACGM)
; The inputs may contain arrays. If both are arrays, then they should be the same size.
;For FORTRAN version, will be easiest to assume that both always contain just
;one element.
; If the location is outside of the model's low-latitude boundary, then a zero
; value is returned, unless the keyword FILL is set, which will be used instead.
; Optional keyword values that are returned are:
;  inside: logical array, set to 1 for each input point within the boundary.
;  outside: logical array, set to 1 for each input point outside of the boundary
;  Bpot= array filled with the magnetic potential values at the same locations.
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc

precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phir

;print,format="('EpotVal: lat=',f8.3,' mlt=',f8.3,' inside=',i3)",$
;  lat,mlt,inside[0]

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)

; print,'EpotVal: ngood=',ngood,' lat=',lat,' mlt=',mlt

  z=fltarr(ngood)
  dobpot=arg_present(bpot)
  if dobpot then z2=fltarr(ngood)

; Dimensions (assuming single lat,mlt): phir[1], phim[1,2]
  phim=phir # replicate(1,maxm) * ((indgen(maxm)+1) ## replicate(1,n_elements(phir)))

; help,phir
; help,phim
; print,format="('phir=',e12.4,' phim=',2e12.4)",phir[0],phim[0,*]

  cospm=cos(phim)
  sinpm=sin(phim)
  For j=0,csize-1 Do Begin
    m=Ms[j]
    if AB[j] then Begin
      Plm=SCPlm(j,colats)

;     print,'Epotval: j=',j,' m=',m,' AB[j]=',AB[j],' colats=',colats,' Plm=',Plm
;     print,format="('Epotval: j=',i3,' plm=',e12.4,' m=',i3)",$
;       j,plm,m

      if m eq 0 then begin
        Z+= Plm*Esphc[j]
        if dobpot then Z2+= Plm*Bsphc[j]
      endif else begin
        mm=m-1
        Z+= Plm*(Esphc[j]*cospm[*,mm]+Esphc[j+1]*sinpm[*,mm])
        if dobpot then Z2+= Plm*(Bsphc[j]*cospm[*,mm]+Bsphc[j+1]*sinpm[*,mm])
        j++
      endelse
    endif
  endfor
; print,'EpotVal: Z=',Z

if ngood EQ nresult Then begin
  if dobpot then bpot=temporary(Z2)
  RETURN,Z
endif
;otherwise:
result=make_array(nresult,VALUE=fill)
if dobpot then begin
  bpot=result
  bpot[goodpts]=Z2
endif
result[goodpts]=Z
Return,result
END
;*********************************************************************************
PRO SetMpFAC,angle,Bt,Tilt,SWVel,SWDen,ALindex,RESET=reset,VERSION=version,YZ=yz,$
  PATH=mypath,SKIPBNDYSET=skipbndyset
;
; Calculate the complete set of SCHA coeficients for the magnetic potential model,
; given an aribitrary IMF angle (degrees from northward toward +Y),
; magnitude Bt (nT), dipole tilt angle (degrees), and
; solar wind velocity (km/sec)
;
;if keyword YZ is set, then By is in place of 'angle' parameter
;and Bz is in place of 'Bt' parameter.
;
COMMON SaveW05scBpot,EX,d1,d2,SCHFits,ALSCHFits
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc
COMMON SetdeltaBscha,cflag,gflag

cflag=(gflag=0) ;clear flags for deltaB calculation; signals that model has changed

if n_elements(SCHFits) eq 0 or keyword_set(Reset) then Begin
  if ~keyword_set(mypath) then mypath=sourcepath('SetMpFAC')
  file=mypath+'W05SCBpot.xdr'
  RESTORE,file
endif

version='W05SC'

if keyword_set(yz) then begin
  default,angle,0.
  default,Bt,-5.
  p1=angle
  p2=bt
  Bt=sqrt(p1^2+p2^2)
  angle=ATAN(p1,p2)*!RADEG
  IF angle LT 0. THEN angle += 360.
endif else begin
  default,angle,180.
  default,Bt,5.
endelse
default,Tilt,0.
default,SWVel,450.
default,SWDen,9.

if ~keyword_set(skipbndyset) Then SetBoundary,angle,Bt,Tilt,SWVel,SWDen,ALindex,PATH=mypath

stilt=sin(tilt*!dtor)
stilt2=stilt^2
sw=Bt*swvel/1000.
swe= (1.-exp(-sw*ex[1]))*sw^ex[0]
C0=1.
SWP=SWvel^2 * SWDen*1.6726E-6
Rang=angle*!DTOR
cosa=cos(Rang)
sina=sin(Rang)
cos2a=cos(2*Rang)
sin2a=sin(2*Rang)
if Bt LT 1. then begin ;remove angle dependency for IMF under 1 nT
  cosa= -1.+bt*(cosa+1.)
  cos2a=1.+bt*(cos2a-1.)
  sina= bt*sina
  sin2a=bt*sin2a
endif

if n_elements(ALindex) GT 0 then Begin
  al=ALindex[0] ;should be a scalar
  cfits=ALSCHFits  ;ALSCHFits=FLTARR(d2,csize)
  A=[C0, swe, stilt, stilt2, swp, al , $
     swe*cosa, stilt*cosa, stilt2*cosa, swp*cosa, al*cosa,$
     swe*sina, stilt*sina, stilt2*sina, swp*sina, al*sina,$
     swe*cos2a,$
     swe*sin2a ]
  d=d2
endif else begin
  cfits=SCHFits  ;;SCHFits=FLTARR(d1,csize)
  ;;A=[C0, swe, swe*cosa, swe*sina, swe*cos2a, swe*sin2a ]
  A=[C0, swe, stilt,stilt2,swp,$
     swe*cosa, stilt*cosa, stilt2*cosa, swp*cosa,  $
     swe*sina, stilt*sina, stilt2*sina, swp*sina,  $
     swe*cos2a,$
     swe*sin2a ]
  d=d1
endelse

Bsphc=fltarr(csize)
for i=0,d-1 do Bsphc += reform(cfits[i,*])*A[i]

if keyword_set(yz) then begin
 ;return inputs to original values
  angle=p1
  Bt=p2
endif

RETURN
END
;*********************************************************************************
FUNCTION BpotVal,lat,mlt,INSIDE=inside,OUTSIDE=outside,FILL_VALUE=fill,EPOT=epot,$
  _REF_EXTRA=extra
;Return the Magnetic Potential at given latitude/MLT,
COMMON SetW05sc,maxM,maxL,Ls,Ms,AB,csize,Esphc,Bsphc

precheckinput,lat,mlt,nresult,ngood,goodpts,inside,colats,phir,_EXTRA=extra
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)
  doepot=arg_present(epot)
  if doepot then z2=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
    m=Ms[j]
    if AB[j] then Begin
      Plm=SCPlm(j,colats)
      if m eq 0 then begin
        Z+= Plm*Bsphc[j]
        if doepot then Z2+= Plm*Esphc[j]
      endif else begin
        mm=m-1
        Z+= Plm*(Bsphc[j]*cospm[*,mm]+Bsphc[j+1]*sinpm[*,mm])
        if doepot then Z2+= Plm*(Esphc[j]*cospm[*,mm]+Esphc[j+1]*sinpm[*,mm])
        j++
      endelse
    endif
  endfor

if ngood EQ nresult Then begin
  if doepot then epot=temporary(Z2)
  RETURN,Z
endif
;otherwise:
result=make_array(nresult,VALUE=fill)
if doepot then begin
  epot=result
  epot[goodpts]=Z2
endif
result[goodpts]=Z
Return,result
END
;********************************************************************
Pro W05 ;calling this causes IDL to compile all routines in this file
Return
End