Алгоритмы палеомагнитных вычислений
ANGLE: Угол между двумя векторами
D1, I1, R=1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
cos(φ) = X1*X2 + Y1*Y2 + Z1*Z2
Входные данные:
- D1 – склонение первого вектора
- I1 – наклонение первого вектора
- D2 – склонение второго вектора
- I2 – наклонение второго вектора
Выходные данные:
- φ – угол между двумя векторами
NORMAL: Нормаль к двум векторам
D1, I1, R=1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
XN = Y1*Z2 – Y2*Z1
YN = Z1*X2 – Z2*X1
ZN = X1*Y2 – X2*Y1
XN, YN, ZN → DN, IN
Входные данные:
- D1 – склонение первого вектора
- I1 – наклонение первого вектора
- D2 – склонение второго вектора
- I2 – наклонение второго вектора
Выходные данные:
- DN – склонение нормали
- IN – наклонение нормали
PROJECT: Проекция одного вектора на другой
D1, I1, R1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
P = R1*(X1*X2 + Y1*Y2 + Z1*Z2)
Входные данные:
- D1 – склонение первого вектора
- I1 – наклонение первого вектора
- R1 – длина первого вектора
- D2 – склонение второго вектора
- I2 – наклонение второго вектора
Выходные данные:
- P – проекция одного вектора на другой
Incl2Lat: Широта из наклонения
tg(Lat) = 0.5*tg(I)
Входные данные:
Выходные данные:
Lat2Incl: Наклонение из широты
tg(I) = 2*tg(Lat)
Входные данные:
Выходные данные:
XYZ2DIR: Сферические координаты из декартовых
p = SQRT(X*X+Y*Y)
R = SQRT(X*X+Y*Y+Z*Z)
D = arccos(X/p)*SIGN(y)
I = arcsin(Z/R)
Входные данные:
- X – координата, направленная на север…
- Y – координата, направленная на восток…
- Z – координата, направленная вниз…
Выходные данные:
- D – склонение…
- I – наклонение…
- R – длина вектора…
DIR2XYZ: Декартовы координаты из сферических
X = R*cos(I)*cos(D)
Y = R*cos(I)*sin(D)
Z = R*sin(I)
Входные данные:
- D – склонение…
- I – наклонение…
- R – длина вектора…
Выходные данные:
- X – координата, направленная на север…
- Y – координата, направленная на восток…
- Z – координата, направленная вниз…
Geo2Strat: Древние координаты из современных
Dg, Ig, R=1 → Xg, Yg, Zg
X= Xg*cos(DipDir)*cos(Dip) +Yg*sin(DipDir)*cos(Dip) +Zg*sin(Dip)
Y= –Xg*sin(DipDir) +Yg*cos(DipDir)
Xs= X*cos(DipDir) –Y*sin(DipDir)
Ys= X*sin(DipDir) +Y*cos(DipDir)
Zs= –Xg*cos(DipDir)*sin(Dip) –Yg*sin(DipDir)*sin(Dip) +Zg*cos(Dip)
Xs, Ys, Zs → Ds, Is
Входные данные:
- Dg – склонение в современных (географических) координатах
- Ig – наклонение в современных (географических) координатах
- DipDir – азимут падения пласта (Dip Direction)
- Dip – угол падения пласта
Выходные данные:
- Ds – склонение в древних (стратиграфических) координатах
- Is – наклонение в древних (стратиграфических) координатах
Strat2Geo: Современные координаты из древних
Ds, Is, R=1 → Xs, Ys, Zs
X= Xs*cos(DipDir)*cos(Dip) +Ys*sin(DipDir)*cos(Dip) –Zs*sin(Dip)
Y= –Xs*sin(DipDir) +Ys*cos(DipDir)
Xg= X*cos(DipDir) –Y*sin(DipDir)
Yg= X*sin(DipDir) +Y*cos(DipDir)
Zg= Xs*cos(DipDir)*sin(Dip) +Ys*sin(DipDir)*sin(Dip) +Zs*cos(Dip)
Xg, Yg, Zg → Dg, Ig
Входные данные:
- Ds – склонение в древних (стратиграфических) координатах
- Is – наклонение в древних (стратиграфических) координатах
- DipDir – азимут падения пласта (Dip Direction)
- Dip – угол падения пласта
Выходные данные:
- Dg – склонение в современных (географических) координатах
- Ig – наклонение в современных (географических) координатах
FindBedding: Определение элементов залегания пласта
DipDir=0;
Dip=0;
if ((Ds<>Dg) or (Is<>Ig)) {
DIR2XYZ(Dg,Ig,1,Xg,Yg,Zg);
DIR2XYZ(Ds,Is,1,Xs,Ys,Zs);
A1=(Ys–Yg)/(Xs-Xg);
B1=Yg-A1*Xg;
DD=-A1-1/A1;
D1=-B1/A1;
D2=B1;
X0=D2/DD;
Y0=D1/DD;
CB=(Xg-X0)*(Xs-X0)+(Yg-Y0)*(Ys-Y0)+Zg*Zs;
CB=CB/SQRT((Xg-X0)*(Xg-X0)+(Yg-Y0)*(Yg-Y0)+Zg*Zg;
CB=CB/SQRT((Xs-X0)*(Xs-X0)+(Ys-Y0)*(Ys-Y0)+Zs*Zs);
DipDir= D360(ATAN(A1)+180);
Dip= ACOS(CB);
Geo2Strat(Xg,Yg,Zg,DirDip,Dip,Xs,Ys,Zs);
XYZ2DIR(Xs,Ys,Zs,Ds_,Is_,Rs_);
if(Angle(Ds,Is,Ds_,Is_)>1e-2) {
DipDir=D360(DipDir-180); } }
Входные данные:
- Ds – склонение в древних (стратиграфических) координатах
- Is – наклонение в древних (стратиграфических) координатах
- Dg – склонение в современных (географических) координатах
- Ig – наклонение в современных (географических) координатах
Выходные данные:
- DipDir – азимут падения пласта
- Dip – угол падения пласта
Dir2Pole: Вычисление координат палеомагнитного полюса
Храмов А.Н. и др. Палеомагнитология. Л.: Недра, 1982. - 312 с. (
стр. 18)
Входные данные:
- D – склонение
- I – наклонение
- SLong – долгота места
- SLat – штрота места
Выходные данные:
- Plong – долгота полюса
- PLat – широта полюса
- PaleoLat – ралеоширота
Pole2Dir: Вычисление направления по известному полюсу
Храмов А.Н. и др. Палеомагнитология. Л.: Недра, 1982. - 312 с. (
стр. 18)
Input data:
- Plong – долгота полюса
- PLat – широта полюса
- SLong – долгота места
- SLat – штрота места
Output data:
- D – склонение
- I – наклонение
- PaleoLat – ралеоширота
PM: Палеомагнитные данные из параметров реконструкции
Incl: из формулы tg(Incl) = 2*tg(PaleoLat)
Decl = D360(Am–Ao)
Pole(Decl, Incl, SLong, SLat, PLong, PLat, Paleolat)
Входные данные:
- PaleoLat – палеоширота
- Am – современный азимут простирания какой-либо структуры
- Ao – азимут простирания той же структуры на реконструкции
- SLong – долгота места
- SLat – широта места
Выходные данные:
- Decl – склонение палеомагнитного вектора
- Incl – наклонение палеомагнитного вектора
- PLong – долгота палеомагнитного полюса
- PLat – широта палеомагнитного полюса
α95: Вычисление α95 из известных N и k
A95=140/SQRT(N*k)
Входные данные:
- k – кучность векторов выборки
- N – объем выборки
Выходные данные:
- A95 – радиус круга доверия
kappa: Вычисление k из известных N и α95
из формулы A95=140/SQRT(N*k)
k=(19600/(A95*A95))/N
Входные данные:
- A95 – радиус круга доверия
- N – объем выборки
Выходные данные:
- k – кучность векторов выборки
N: Вычисление N из известных k и α95
из формулы A95=140/SQRT(N*k)
N=(19600/(A95*A95))/k
Входные данные:
- A95 – радиус круга доверия
- k – кучность векторов выборки
Выходные данные:
R(r): Вектор-результант из кучности и объема выборки
из формулы k = (N–1)/(N–R)
R = N – (N – 1)/k
r = R/N
Входные данные:
- k – кучность векторов выборки
- N – объем выборки
Выходные данные:
- R – вектор-результант
- r – нормированный вектор-результант
Sum: Сумма двух векторов
N=N1+N2
D1, I1, R1 → X1, Y1, Z1
D2, I2, R2 → X2, Y2, Z2
X=X1+X2
Y=Y1+Y2
Z=Z1+Z2
X, Y, Z → Dm, Im, R
k=(N-1)/(N-R)
a95=140/SQRT(k*N)
Входные данные:
- D1 – склонение первого вектора
- I1 – наклонение первого вектора
- R1 – длина первого вектора
- N1 – объем первой выборки
- D2 – склонение второго вектора
- I2 – наклонение второго вектора
- R2 – длина второго вектора
- N2 – объем второй выборки
Выходные данные:
- Dm – склонение суммарного вектора
- Im – наклонение сумарного вектора
- R – длина суммарного вектора
- N – кучность
- k – радиус круга доверия на уровне 95%
- A95 – объем суммарной выборки
Примечание:
- R1, N1, R2 и N2 – могут быть равны и единице :)
Euler: Вращение на сфере вокруг полюса Эйлера
Корн Г., Корн Т. Справочник по математике (ф-лы 14.10-6 – 14.10-8)
Входные данные:
- L1 – долгота полюса до поворота
- F1 – широта полюса до поворота
- LE – долгота полюса Эйлера
- FE – широта полюса Эйлера
- FI – угол поворота
Выходные данные:
- L2 – долгота полюса после поворота
- F2 – широта полюса после поворота
RevTest: Тест обращения (McFadden, McElhinny, 1990)
McFadden P.L., McElhinny M.W. Classification of the reversal test in palaeomagnetism. Geophys. J. Int. (1990) 103, 725-729
Входные данные:
- D1 – склонение среднего первой выборки
- I1 – наклонение среднего первой выборки
- k1 – кучность первой выборки
- N1 – объем первой выборки
- D1 – склонение среднего второй выборки
- I1 – наклонение среднего второй выборки
- k1 – кучность второй выборки
- N1 – объем второй выборки
Выходные данные:
- Ratio – отношение кучностей k1/k2
- Rc – критическое значение для отношения кучностей
- Dm – склонение суммарного вектора
- Im – наклонение суммарного вектора
- R – длина суммарного вектора
- k – кучность суммарной выборки
- a95 – радиус круга доверия на уровне 95%
- N – объем суммарной выборки
- Gamma – угол между средними векторами (гамма)
- Gc – критическое значение для угла
Algorithms for paleomagnetic calculation
Angle: Angle between two vectors
D1, I1, R=1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
cos(φ) = X1*X2 + Y1*Y2 + Z1*Z2
Input data:
- D1 – declination of first vector
- I1 – inclination of first vector
- D2 – declination of second vector
- I2 – inclination of second vector
Output data :
- φ – angle between two vectors
Normal: Perpendicular (normal) vector to two vectors
D1, I1, R=1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
Xn = Y1*Z2 – Y2*Z1
Yn = Z1*X2 – Z2*X1
Zn = X1*Y2 – X2*Y1
Xn, Yn, Zn → Dn, In
Input data
- D1 – declination of first vector
- I1 – inclination of first vector
- D2 – declination of second vector
- I2 – inclination of second vector
Output data
- Dn – declination of normal vector
- In – inclination of normal vector
Project: Projection of a vector to the given vector
D1, I1, R1 → X1, Y1, Z1
D2, I2, R=1 → X2, Y2, Z2
p = R1*(X1*X2 + Y1*Y2 + Z1*Z2)
Input data:
- D1 – declination of first vector
- I1 – inclination of first vector
- R1 – length of first vector
- D2 – declination of second vector
- I2 – inclination of second vector
Output data:
- p – projection of a vector to the given vector
Incl2Lat: Latitude from inclination
tg(Lat) = 0.5*tg(I)
Input data:
- I – inclination of vector
Output data:
Lat2Incl: Inclination from latitude
tg(I) = 2*tg(Lat)
Input data:
Output data:
- I – inclination of vector
XYZ2DIR: Spherical coordinates from Cartesian
p = SQRT(X*X+Y*Y)
R = SQRT(X*X+Y*Y+Z*Z)
D = arccos(X/p)*SIGN(y)
I = arcsin(Z/R)
Input data:
- X – north Cartesian coordinate
- Y – east Cartesian coordinate
- Z – down Cartesian coordinate
Output data:
- D – declination
- I – inclination
- R – verctor length
DIR2XYZ: Cartesian coordinates from Spherical
X = R*cos(I)*cos(D)
Y = R*cos(I)*sin(D)
Z = R*sin(I)
Input data:
- D – declination
- I – inclination
- R – verctor length
Output data:
- X – north Cartesian coordinate
- Y – east Cartesian coordinate
- Z – down Cartesian coordinate
Geo2Strat: Stratigraphic coordinates from geographic
Dg, Ig, R=1 → Xg, Yg, Zg
X= Xg*cos(DipDir)*cos(Dip) +Yg*sin(DipDir)*cos(Dip) +Zg*sin(Dip)
Y= –Xg*sin(DipDir) +Yg*cos(DipDir)
Xs= X*cos(DipDir) –Y*sin(DipDir)
Ys= X*sin(DipDir) +Y*cos(DipDir)
Zs= –Xg*cos(DipDir)*sin(Dip) –Yg*sin(DipDir)*sin(Dip) +Zg*cos(Dip)
Xs, Ys, Zs → Ds, Is
Input data:
- Dg – declination in geographic coordinates
- Ig – inclination in geographic coordinates
- DipDir – dip direction
- Dip – dip
Output data:
- Ds – declination in stratigraphic coordinates
- Is – inclination in stratigraphic coordinates
Strat2Geo: Geographic coordinates from stratigraphic
Ds, Is, R=1 → Xs, Ys, Zs
X= Xs*cos(DipDir)*cos(Dip) +Ys*sin(DipDir)*cos(Dip) –Zs*sin(Dip)
Y= –Xs*sin(DipDir) +Ys*cos(DipDir)
Xg= X*cos(DipDir) –Y*sin(DipDir)
Yg= X*sin(DipDir) +Y*cos(DipDir)
Zg= Xs*cos(DipDir)*sin(Dip) +Ys*sin(DipDir)*sin(Dip) +Zs*cos(Dip)
Xg, Yg, Zg → Dg, Ig
Input data:
- Ds – declination of vector in stratigraphic coordinates
- Is – inclination of vector in stratigraphic coordinates
- DipDir – dip direction
- Dip – dip
Output data:
- Dg – declination of vector in geographic coordinates
- Ig – inclination of vector in geographic coordinates
FindBedding: Estimation of bedding geometries
DipDir=0;
Dip=0;
if ((Ds<>Dg) or (Is<>Ig)) {
DIR2XYZ(Dg,Ig,1,Xg,Yg,Zg);
DIR2XYZ(Ds,Is,1,Xs,Ys,Zs);
A1=(Ys–Yg)/(Xs-Xg);
B1=Yg-A1*Xg;
DD=-A1-1/A1;
D1=-B1/A1;
D2=B1;
X0=D2/DD;
Y0=D1/DD;
CB=(Xg-X0)*(Xs-X0)+(Yg-Y0)*(Ys-Y0)+Zg*Zs;
CB=CB/SQRT((Xg-X0)*(Xg-X0)+(Yg-Y0)*(Yg-Y0)+Zg*Zg;
CB=CB/SQRT((Xs-X0)*(Xs-X0)+(Ys-Y0)*(Ys-Y0)+Zs*Zs);
DipDir= D360(ATAN(A1)+180);
Dip= ACOS(CB);
Geo2Strat(Xg,Yg,Zg,DirDip,Dip,Xs,Ys,Zs);
XYZ2DIR(Xs,Ys,Zs,Ds_,Is_,Rs_);
if(Angle(Ds,Is,Ds_,Is_)>1e-2) {
DipDir=D360(DipDir-180); } }
Input data:
- Ds – declination of vector in stratigraphic coordinates
- Is – inclination of vector in stratigraphic coordinates
- Dg – declination of vector in geographic coordinates
- Ig – declination of vector in geographic coordinates
Output data:
- DipDir – dip direction
- Dip – dip
Dir2Pole: Calculating of virtual geomagnetic pole
Pole2Dir: Calculating a direction from virtual geomagnetic pole
PM: Paleomagnetic data from reconstruction parameters
Incl: from formula tg(Incl) = 2*tg(PaleoLat)
Decl = D360(Am–Ao)
Pole(Decl, Incl, SLong, SLat, PLong, PLat, Paleolat)
Input data:
- PaleoLat – paleolatitude
- Am – modern strike azimuth of the object
- Ao – ancient strike azimuth of the same object
- SLong – site longitude
- SLat – site latitude
Output data:
- Decl – declination of paleomagnetic vector
- Incl – inclination of paleomagnetic vector
- PLong – pole longitude
- PLat – pole latitude
α95: Calculating of circular confidence limit from k and N
A95=140/SQRT(N*k)
Input data:
- k – precision parameter (kappa)
- N – number of vectors
Output data:
- A95 – circular confidence limit
kappa: Calculating of precision parameter from N and α95
from formulas A95=140/SQRT(N*k)
k=(19600/(A95*A95))/N
Input data:
- A95 – circular confidence limit
- N – number of vectors
Output data:
- k – precision parameter (kappa)
N: Calculating of number of vectors from k and α95
from formulas A95=140/SQRT(N*k)
N=(19600/(A95*A95))/k
Input data:
- A95 – circular confidence limit
- k – precision parameter (kappa)
Output data:
R(r): Vector-resultant from from k (kappa) and N
from formulas k = (N–1)/(N–R)
R = N – (N – 1)/k
r = R/N
Input data:
- k – precision parameter (kappa)
- N – number of vectors
Output data:
- R – normalized vector-resultant
- r – normalized vector-resultant
Sum: Sum of two vectors
N=N1+N2
D1, I1, R1 → X1, Y1, Z1
D2, I2, R2 → X2, Y2, Z2
X=X1+X2
Y=Y1+Y2
Z=Z1+Z2
X, Y, Z → Dm, Im, R
k=(N-1)/(N-R)
A95=140/SQRT(k*N)
Input data:
- D1 – declination of first vector
- I1 – inclination of first vector
- R1 – length of first vector
- N1 – size of first sample
- D2 – declination of second vector
- I2 – inclination of second vector
- R2 – length of second vector
- N2 – size of second sample
Output data:
- Dm – declination of total vector
- Im – inclination of total vector
- R – length of total vector
- N – size of total sample
- k – precision parameter (kappa)
- A95 – circular confidence limit
Note:
- R1, N1, R2 & N2 – may be equal to unity :)
Euler: Rotation on the sphere around the Euler pole
Mathematical Handbook for Scientists and Engineers. by Granino A. Korn, Theresa M. Korn (formulas 14.10-6 – 14.10-8)
Input data:
- L1 – longitude of pole before rotation
- F1 – latitude of pole before rotation
- LE – longitude of Euler pole
- FE – latitude of Euler pole
- FI – angle of rotation
Output data:
- L2 – longitude of pole after rotation
- F2 – latitude of pole after rotation
RevTest: Reversal test (McFadden, McElhinny, 1990)
McFadden P.L., McElhinny M.W. Classification of the reversal test in palaeomagnetism. Geophys. J. Int. (1990) 103, 725-729
Input data:
- D1 – total declination of first sample
- I1 – total inclination of first sample
- k1 – precision parameter (kappa) of first sample
- N1 – number of vectors of first sample
- D1 – total declination of second sample
- I1 – total inclination of second sample
- k1 – precision parameter (kappa) of second sample
- N1 – number of vectors of second sample
Output data:
- Ratio – ratio k1/k2
- Rc – critical value for this ratio
- Dm – declination of total vector
- Im – inclination of total vector
- R – length of total vector
- k – precision parameter (kappa) of total sample
- a95 – circular confidence limit
- N – size of total sample
- Gamma – angle beetween mean vectors of two samples
- Gc – critical value for this angle