Перевод географических координат в геомагнитные...

yurimedvedev

Как делать?
Есть ли однозначное соответствие между географической широтой/долготой и склонением/наклонением? А между склонением/наклонением и магнитной широтой/долготой? Я запутался. Второй вопрос интересует больше всего.
Вопрос родился откуда? Искал модель геомагнитного поля (программу) IGRF, нашел несколько разных, которые считают разные параметры, хочу объединить в одну и проверить насколько они правильно работают.

yurimedvedev

Вроде, нашел функцию, которая это делает.
Знатоки, поясните, что же хотел сказать аффтар тут:
http://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/iri90/iri90_storm/irifu12.for
 
C************************************************************
C*************** EARTH MAGNETIC FIELD ***********************
C**************************************************************
C
C
C
SUBROUTINE GGM(ART,LONG,LATI,MLONG,MLAT)
C CALCULATES GEOMAGNETIC LONGITUDE (MLONG) AND LATITUDE (MLAT)
C FROM GEOGRAFIC LONGITUDE (LONG) AND LATITUDE (LATI) FOR ART=0
C AND REVERSE FOR ART=1. ALL ANGLES IN DEGREE.
C LATITUDE:-90 TO 90. LONGITUDE:0 TO 360 EAST.
INTEGER ART
REAL MLONG,MLAT,LONG,LATI
COMMON/CONST/FAKTOR
ZPI=FAKTOR*360.
CBG=11.4*FAKTOR
CI=COS(CBG)
SI=SIN(CBG)
IF(ART.EQ.0) GOTO 10
CBM=COS(MLAT*FAKTOR)
SBM=SIN(MLAT*FAKTOR)
CLM=COS(MLONG*FAKTOR)
SLM=SIN(MLONG*FAKTOR)
SBG=SBM*CI-CBM*CLM*SI
IF(ABS(SBG).GT.1.) SBG=SIGN(1.,SBG)
LATI=ASIN(SBG)
CBG=COS(LATI)
SLG=(CBM*SLM)/CBG
CLG=(SBM*SI+CBM*CLM*CI)/CBG
IF(ABS(CLG).GT.1.) CLG=SIGN(1.,CLG)
LONG=ACOS(CLG)
IF(SLG.LT.0.0) LONG=ZPI-LONG
LATI=LATI/FAKTOR
LONG=LONG/FAKTOR
LONG=LONG-69.8
IF(LONG.LT.0.0) LONG=LONG+360.0
RETURN
10 YLG=LONG+69.8
CBG=COS(LATI*FAKTOR)
SBG=SIN(LATI*FAKTOR)
CLG=COS(YLG*FAKTOR)
SLG=SIN(YLG*FAKTOR)
SBM=SBG*CI+CBG*CLG*SI
IF(ABS(SBM).GT.1.) SBM=SIGN(1.,SBM)
MLAT=ASIN(SBM)
CBM=COS(MLAT)
SLM=(CBG*SLG)/CBM
CLM=(-SBG*SI+CBG*CLG*CI)/CBM
IF(ABS(CLM).GT.1.) CLM=SIGN(1.,CLM)
MLONG=ACOS(CLM)
IF(SLM.LT..0) MLONG=ZPI-MLONG
MLAT=MLAT/FAKTOR
MLONG=MLONG/FAKTOR
RETURN
END

yurimedvedev

Тут есть два параметра, 11.4 градуса и 69.8. Они наверняка связаны с положением магнитного полюса относительно географического. Однако они постоянно меняются и это не есть хорошо. Нужны эти параметры в зависимости от времени.

meles

Глюк, попробуй поюзать MapInfo
может быть там есть
(непонятно, тебе нужен алгоритм или результат)

yurimedvedev

Попробую.
Мне нужен алгоритм. На любом языке, хотя фортран уже стал предпочительным.

k1a2r3t4a

SUBROUTINE TRANS (UT,IDAY,iyear,tpsi,BD)
***************************************************************
* Calculation of the geomagnetic dipole tilt angle and .
* matrices of transition from the geographic coordinates .
* (Cartesian) to the GSM-coordinates (G2GSM(3,3) in the .
* COMMON BLOCK /TRAN/). .
* _ _ .
* Xgsm = G2GSM*Xgeogr .
* .
* ALPHA1 is the angle between the Earth's axis and the .
* dipole moment; .
* ALPHA2 is the angle between the Earth's axis and .
* the normal to the ecliptic plane; .
* PHIM is the angle between the midnight meridian plane .
* and the meridional plane; .
* PHISE is the angle between the Earth-Sun line and .
* the projection of the Earth's axis on the .
* ecliptic plane; .
* PSI is the tilt angle of the geomagnetic dipole; .
* B is the Sun's deflection; .
* UT is universal time; .
* IDAY is the day number; .
* B1 is western longitude of the noon meridian; .
* B2 is the angle between the noon geogr. meridian .
* and the Y=0 plane in the GSM-coordinates. .
* Written by I. Alexeev .
* .
* Acknowledgment: .
* Program SUN written by Gilbert D. Mead is used .
* to determine some parameters dependent on the Earth .
* and Sun mutual position (SLONG, sind, and cosd from .
* COMMON BLOCK /ddd/). .
***************************************************************
COMMON/TRAN/g2gsm(3,3)
common /ddd/ sind,cosd
dimension gauss(3)
Ihour=int(ut)
dmin=(ut-ihour)*60.
min=int(dmin)
isec=intdmin-min)*60.)
call SUN(IYeaR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
P= 3.1415926/180.
pi2=3.1415926/2
call dipgarm(iyear, gauss)
G1=GAUSS(2)
H1=GAUSS(3)
PD=G1**2+H1**2
G10=GAUSS(1)
BD=-SQRT(G10*G10+PD)
ALPHA1= ATAN(-1.*SQRT(PD)/G10)
PHINP= ATAN(H1/G1)
ALPHA2= 23.4419*P
PHIM= P*(UT*15+phinp/p)
PHISE= pi2-slong+9.924E-5

B1=P*15.*(UT-12)
SB= SIN(ALPHA2)*COS(PHISE)
CB= SQRT(1-SB*SB)
SB1=SIN(B1)
CB1=COS(B1)
SPSI= -SB*COS(ALPHA1) + CB*SIN(ALPHA1)*COS(PHIM)
CPSI= SQRT(1-SPSI*SPSI)
psi=asin(spsi)/p
CB2=(COS(ALPHA1)+SB*SPSI)/CPSI/CB
SB2=SQRT(abs(1-CB2*CB2
IF(PHIM.LE.0..OR.PHIM.GE.3.1415926) SB2=-SB2
tpsi=psi
g2gsm(1,1)=cb1*cb
g2gsm(1,2)=-sb1*cb
g2gsm(1,3)=sb
g2gsm(2,1)=sb1*cb2-cb1*sb*sb2
g2gsm(2,2)=cb1*cb2+sb1*sb*sb2
g2gsm(2,3)=cb*sb2
g2gsm(3,1)=-sb1*sb2-cb1*sb*cb2
g2gsm(3,2)=-cb1*sb2+sb1*sb*cb2
g2gsm(3,3)=cb*cb2
RETURN
END

FUNCTION IDD(iy,mo,id)
********************************************************************
* Calculation of the day number in a year *
* INPUT PARAMETERS: year (IY month (MO day in the month (ID) *
* Written by V. Kalegaev *
********************************************************************
II=0
IF (MO.EQ.1) GOTO 4
5 DO 10 M=1,MO-1
GOTO (1,2,1,3,1,3,1,1,3,1,3) M
1 II=II+31
GOTO 10
2 II=II+28
GOTO 10
3 II=II+30
10 CONTINUE
4 II=II+ID
if (mod(iy,100).eq.0.and.mod(iy,400).ne.0) goto 6
IF (MOD(IY,4).EQ.0.AND.II.GT.59.AND.MO.GT.2) II=II+1
6 IDD=II
RETURN
END
2. aplha1 = 11,43 град - угол наклона геомагнитного диполя к оси вращения Земли
-69б16 град - географическая широта северного геомагнитного полюса.
psi - угол наклона земного диполя к плоскости перпендикулярной лини Земля- Солнце ( >0 - зима в северном полушарии, <0 - лето). Psi зависит от времени суток и времени года.
sin(psi) = - cos(alpha1)sin(betta2) + sin(aplha1)cos(betta2)cos(phi_m
где betta2 -склонение Солнца, зависящее от номенра дня в году IDAY (в течении суток изменяется мало).
phi_m - угол между плоскостью полуночного меридиана и меридиональной плоскостью, содержащей северный магнитный полюс, зависит от м ирового времени UT:
phi_m= pi*(UT/12-69,16/180)
В течение года psi меняется в пределах +_35 град, а в течение суток - в пределах +_11б43 град относитльено некотрого значения, равного - betta2, зависящего от IDAY
Оставить комментарий
Имя или ник:
Комментарий: