PRO c_r_to_a, timearray, timereal

  dayofmon = [31,28,31,30,31,30,31,31,30,31,30,31]

  timearray = intarr(6)

  speryear = double(31536000.0)
  sperday  = double(86400.0)
  sperhour = double(3600.0)
  spermin  = double(60.0)

  numofyears = floor(timereal/speryear)
  if (numofyears+65) mod 4 eq 0 then dayofmon(1) = dayofmon(1) + 1
  numofdays = floor((timereal mod speryear)/sperday)
  numofleap = floor(numofyears / 4)
  numofdays = numofdays - numofleap
  if numofdays lt 0 then begin
    if (numofyears+65) mod 4 eq 0 then dayofmon(1) = dayofmon(1) - 1
    numofyears = numofyears - 1
    numofdays = numofdays + numofleap + 365
    if (numofyears+65) mod 4 eq 0 then dayofmon(1) = dayofmon(1) + 1
    numofleap = floor(numofyears / 4)
    numofdays = numofdays - numofleap
  endif
  numofhours = floor((timereal mod sperday)/sperhour)
  numofmin = floor((timereal mod sperhour)/spermin)
  numofsec = floor(timereal mod spermin)

  numofmon = 0

  while numofdays ge dayofmon(numofmon) do begin

    numofdays = numofdays - dayofmon(numofmon)
    numofmon = numofmon + 1

  endwhile

  timearray(0) = numofyears + 1965
  timearray(1) = numofmon + 1
  timearray(2) = numofdays + 1
  timearray(3) = numofhours
  timearray(4) = numofmin
  timearray(5) = numofsec

  RETURN

END

filelist = findfile("b*.save")

filename = filelist(0)

print, "restoring file : ", filename

restore, filename

nFields = n_elements(fields)

for i=0,nFields-1 do print, string(i,format="(i2)")+". ",fields(i)

iTimes = fix(randomu(s,5)*1440)

nLats = n_elements(data(0,0,0,*))
nMLTs = n_elements(data(0,0,*,0))

iLats = fix(randomu(s,5)*nLats)
iMLTs = fix(randomu(s,5)*nMLTs)

iEEast_ = 7
iENorth_ = 8

for i = 0, 4 do begin

    c_r_to_a, iArray, time(iTimes(i))

    for iLoc = 0, 4 do begin

        print, "East/Northward E-field at ",iArray
        print, "         and  location ",lats(iLats(iLoc)),$
          mlts(iMLTs(iLoc))," : "
        print, $
          data(iTimes(i),iEEast_, iMLTs(iLoc), iLats(iLoc)), " V/m", $
          data(iTimes(i),iENorth_, iMLTs(iLoc), iLats(iLoc)), " V/m"

    endfor
endfor

print, "Computing Field-align currents"

iEField_ = -1
iPed_    = -1

for iField = 0, nFields-1 do begin
    fields(iField) = mklower(fields(iField))
    if (strpos(fields(iField),'electric field') gt -1 and $
        iEField_ eq -1) then iEField_ = iField
    if (strpos(fields(iField),'ped') gt -1 and $
        strpos(fields(iField),'total') gt -1 and $
        iPed_ eq -1) then iPed_ = iField
endfor


nTimes = n_elements(time)

fac_all = fltarr(nTimes, nMlts, nLats)

    theta = (90.0 - lats)*!dtor
    phi   = mlts * !pi/12.0
    re    = 6372000.0

    fac = fltarr(nMlts, nLats)
    for iTime = 0, nTimes-1 do begin
        ee = reform(data(iTime, iEField_+0, *, *))
        en = reform(data(iTime, iEField_+1, *, *))
        sp = reform(data(iTime, iPed_,      *, *))
        je = ee * sp
        jn = -en * sp 
        for iMlts=0,nMlts-1 do jn(iMlts,*) = jn(iMlts,*) * sin(theta)
        fac(*,*) = 0.0

        ; North Component
        for iMlts=0,nMlts-2 do begin
            fac(iMlts,1:nLats-2) = $
              (jn(iMlts,2:nLats-1) - jn(iMlts,0:nLats-3)) / $
              (theta(2:nLats-1)-theta(0:nLats-3))
            fac(iMlts,nLats-1)   = $
              (jn(iMlts,nLats-1) - jn(iMlts,nLats-2)) / $
              (theta(nLats-1)-theta(nLats-2))
        endfor

        ; East Component
        for iLat = 1, nLats-1 do begin
            fac(1:nMlts-2, iLat) = fac(1:nMlts-1, iLat) + $
              (je (2:nMlts-1,iLat) - je (0:nMlts-3,iLat)) / $
              (phi(2:nMlts-1)      - phi(0:nMlts-3))
            ; across the 0-24 border
            fac(0, iLat) = fac(0, iLat) + $
              (je (1,iLat) - je (nMlts-2,iLat)) / $
              (phi(2)      - phi(0))
            fac(*, iLat) = fac(*, iLat)/(re*sin(theta(iLat)))
        endfor

        fac(nMlts-1, *) = fac(0, *)
        fac_all(iTIme, *, *) = -fac*1.0e6

    endfor



end