pro geo_to_apex, geoglat, geoglon, apexlat, apexlon, apexfile ; Transform geographic coordinates to apex magnetic coordinates. ; This requires a data file (supplied by Art Richmond) which contains ; the apex coordinates for all geographic coordinates on a 1x1 degree ; grid. The routine works for all latitudes and longitudes, and the ; data file is valid only for the epoch 1997. ; keep the data in a common block, i.e. read the file only once common geo_to_apex_data, data if n_elements(data) le 0 then begin alon = fltarr(361,181) alat = fltarr(361,181) openr, unit, apexfile, /get_lun, /xdr readu, unit, alat readu, unit, alon free_lun, unit endif ; the interpolation requires longitude 0 to 360 ; Map the line segment (-180,180) into a circle in the ; complex plane, perform the interpolation, map back to ; the original line segment salon = sin(!DTOR*alon) calon = cos(!DTOR*alon) sapexlon = interpolate(salon, ((geoglon+360) mod 360), geoglat + 90) capexlon = interpolate(calon, ((geoglon+360) mod 360), geoglat + 90) apexlon = atan(sapexlon,capexlon)/!DTOR apexlat = interpolate(alat, ((geoglon+360) mod 360), geoglat + 90) ; the interpolation requires longitude 0 to 360 ;apexlat = interpolate(alat, ((geoglon+360) mod 360), geoglat + 90) ;apexlon = interpolate(alon, ((geoglon+360) mod 360), geoglat + 90) ; note: interpolation causes problems at the alon=180 border.... ;apexlat = alat(((geoglon+360) mod 360), geoglat + 90) ;apexlon = alon(((geoglon+360) mod 360), geoglat + 90) end