地理座標から地磁気座標への変換

Question

地理座標と地磁気座標の間の変換を試みています。私は次のPrologスクリプトを見つけましたが、私はそれを自分で変換するのに十分理解していません。ターゲット言語はJavaですが、理解できるものは何でもいいです（C、Python、VBなど）。

http://idlastro.gsfc.nasa.gov/ftp/pro/astro/geo2mag.pro

誰かがこのスクリプトの変換を助けるか、それは（それらの配列操作が私に不可解されている）やっているまさに説明できるいずれかの場合、私はそれを本当に感謝します。

おかげ

Answer 1

私はPythonコードにそれを作って、このウェブサイトhttp://wdc.kugi.kyoto-u.ac.jp/igrf/gggm/index.htmlで検証してみました。私は

1. 磁極が、私は1995年の値を使用するために、上記の計算を設定しても年1995年
2. のために、私はかなり正しい一致を取得しないということであることがわかりました。

京都、日本（35N、135.45W）の値を使用しました。 Webページが計算されました（25.18、-155.80）。私は得た（25.33580652、-155.82724011）。私は、アプリケーションに応じて...

import numpy as np 

from numpy import pi, cos, sin, arctan2, sqrt, dot 
def geo2mag(incoord): 
    """geographic coordinate to magnetic coordinate: 

     incoord is numpy array of shape (2,*) 
     array([[glat0,glat1,glat2,...], 
      [glon0,glon1,glon2,...]) 
     where glat, glon are geographic latitude and longitude 
     (or if you have only one point it is [[glat,glon]]) 

     returns 
     array([mlat0,mlat1,...], 
      [mlon0,mlon1,...]]) 
     """ 

    # SOME 'constants'... 
    lon = 288.59 # or 71.41W 
    lat = 79.3 
    r = 1.0 

    # convert first to radians 
    lon, lat = [x*pi/180 for x in lon,lat] 

    glat = incoord[0] * pi/180.0 
    glon = incoord[1] * pi/180.0 
    galt = glat * 0. + r 

    coord = np.vstack([glat,glon,galt]) 

    # convert to rectangular coordinates 
    x = coord[2]*cos(coord[0])*cos(coord[1]) 
    y = coord[2]*cos(coord[0])*sin(coord[1]) 
    z = coord[2]*sin(coord[0]) 
    xyz = np.vstack((x,y,z)) 

    # computer 1st rotation matrix: 
    geo2maglon = np.zeros((3,3), dtype='float64') 
    geo2maglon[0,0] = cos(lon) 
    geo2maglon[0,1] = sin(lon) 
    geo2maglon[1,0] = -sin(lon) 
    geo2maglon[1,1] = cos(lon) 
    geo2maglon[2,2] = 1. 
    out = dot(geo2maglon , xyz) 

    tomaglat = np.zeros((3,3), dtype='float64') 
    tomaglat[0,0] = cos(.5*pi-lat) 
    tomaglat[0,2] = -sin(.5*pi-lat) 
    tomaglat[2,0] = sin(.5*pi-lat) 
    tomaglat[2,2] = cos(.5*pi-lat) 
    tomaglat[1,1] = 1. 
    out = dot(tomaglat , out) 

    mlat = arctan2(out[2], 
      sqrt(out[0]*out[0] + out[1]*out[1])) 
    mlat = mlat * 180/pi 
    mlon = arctan2(out[1], out[0]) 
    mlon = mlon * 180/pi 

    outcoord = np.vstack((mlat, mlon)) 
    return outcoord 

if __name__ == '__main__': 
    mag = geo2mag(np.array([[79.3,288.59]]).T).T 
    print mag # should be [90,0] 

    mag = geo2mag(np.array([[90,0]]).T).T 
    print mag # should be [79.3,*] 

    mag = geo2mag(np.array([ 
     [79.3,288.59], 
     [90,0] 
     ]).T).T 

    print mag # should be [ [90,0]. [79.3,*] ] 

    # kyoto, japan 
    mag = geo2mag(np.array([[35.,135.45]]).T).T 
    print mag # should be [25.18, -155.80], according to 
       # this site using value for 1995 
       # http://wdc.kugi.kyoto-u.ac.jp/igrf/gggm/index.html 

Answer 2

これは、実際の利用可能ならば、完全にわからない地磁気座標は地球の双極子磁場のマッピングされているので、だから、高度はこの座標変換で重要な変数であることができます。

Pythonでは、SpacePy（http://sourceforge.net/projects/spacepy/）を使用して地理座標を地磁気（またはその逆）に簡単に変換できます。

は、Javaに変換するためのソースコードを探しているので、SpacePyは場合には、Pythonでは提供されていますソースがFortranの国際放射線帯環境モデリング（IRBEM）ライブラリー、（http://irbem.svn.sourceforge.net/viewvc/irbem/web/index.html）

を実施しています他の人は簡単な解決策を探しています：

import spacepy.coordinates as coord 
from spacepy.time import Ticktock 
import numpy as np 
def geotomag(alt,lat,lon): 
    #call with altitude in kilometers and lat/lon in degrees 
    Re=6371.0 #mean Earth radius in kilometers 
    #setup the geographic coordinate object with altitude in earth radii 
    cvals = coord.Coords([np.float((alt/Re+Re))/Re,np.float(lat),np.float(lon)], 'GEO', 'sph',['Re','deg','deg']) 
    #set time epoch for coordinates: 
    cvals.ticks=Ticktock(['2012-01-01T12:00:00'], 'ISO') 
    #return the magnetic coords in the same units as the geographic: 
    return cvals.convert('MAG','sph')