c############################################################################ c c A R A L E V - Main program for the Arason-Levi Maximum Likelihood c Estimates (MLE) for Inclination-only Data. c The method is described in a report and at: c Arason, P. and S. Levi (2006), The Maximum Likelihood Solution for c Inclination-only Data, Eos Trans. AGU, 87(52), Fall Meet. Suppl., c Abstract GP21B-1312. c See also: http://www.vedur.is/~arason/paleomag c c Author: c Dr. Pordur Arason, arason@vedur.is c Icelandic Meteorological Office, Bustadavegi 9, c IS-150 Reykjavik, Iceland c c History: c Winter 1991 Vision that this could be done c Summer 2004 Decision to implement c July 2006 Implementation and testing c Nov 2006 Comment clarifications c Apr 2009 Modification of difficult situations c August 2009 Comment clarifications c c############################################################################ program ARALEV_PR real xinc(10000) call OPT np = 0 dowhile ( .TRUE. ) np = np + 1 read ( *, '(F20.0)', ERR=29, END=19 ) xinc(np) enddo 19 np = np - 1 call AVERAG ( xinc, np, ainc, ak, t63, a95 ) write ( *, '(I6,F8.3,G14.6,2F8.3)' ) np, ainc, ak, t63, a95 call ARALEV ( xinc, np, ainc, ak, t63, a95, ierr ) write ( *, '(I6,F8.3,G14.6,2F8.3,I3)' ) np, ainc, ak, t63, a95, ierr stop 29 write ( 0, '('' ERROR in reading data'')' ) stop end c############################################################################ c c O P T - Runtime options c c############################################################################ subroutine OPT real*4 para character oa*50 call A_NARG ( na ) io = 0 do i=1, na io = io + 1 if ( io .gt. na ) goto 99 call A_GETARG ( io, oa ) if ( oa(1:2) .eq. '-h' ) then write ( 0, '( 1 ''aralev - A program to calculate Arason-Levi'', 1 '' statistics.''/ 1 ''Author: Žóršur Arason - Vešurstofunni - November 2006''/ 1 ''Usage: aralev < INCL-list > stdout''/ 1 ''Output: N Incl Kappa Theta-63 Alpha-95 ''/ 1 ''Options: ''/ 1 '' -h : Help'')') stop endif enddo 99 return end c############################################################################ c c A R A L E V - The Arason-Levi Maximum Likelihood Estimates (MLE) c for Inclination-only Data. c The method is described in a report and at: c Arason, P. and S. Levi (2006), The Maximum Likelihood Solution for c Inclination-only Data, Eos Trans. AGU, 87(52), Fall Meet. Suppl., c Abstract GP21B-1312. c See also: http://www.vedur.is/~arason/paleomag c c xinc input Vector of inclinations (in degrees) c n input Number of data in the sample c ainc output MLE of mean inclination (in degrees) c ak output MLE of precision parameter (kappa) c t63 output Angular standard deviation (Theta-63) c a95 output 95% confidence limits of the mean (Alpha-95) c ierr output Error indicator: c 0: no problem c 1: convergence problem of selected solution c 2: failed robustness check c 3: convergence and robustness problems c c Author: c Dr. Pordur Arason, arason@vedur.is c Icelandic Meteorological Office, Bustadavegi 9, c IS-150 Reykjavik, Iceland c c History: c Winter 1991 Vision that this could be done c Summer 2004 Decision to implement c July 2006 Implementation and testing c Nov 2006 Comment clarifications c Apr 2009 Modification of difficult situations c August 2009 Comment clarifications c c############################################################################ subroutine ARALEV ( xinc, n, ainc, ak, t63, a95, ierr ) real xinc(10000), ainc, ak, t63, a95 integer n integer ierr, ie1, ie2, ie3 real*8 th(10000), dr, dn, t63max, a95max real*8 s, s2, c, x, rt, rk, rt1, rk1, co, dk, dt real*8 the1, the2, the3, the4, akap1, akap2, akap3, akap4 real*8 xl, xl1, xl2, xl3, xl4 real*8 AL1, AL2, COTH, XLIK c Set constants dr = 0.0174532925199433D0 ! Degrees to radians (pi/180) t63max = 105.070062145D0 ! 63 % of a sphere. a95max = 154.158067237D0 ! 95 % of a sphere. dn = n ierr = 1 c Check for illegal use if ( n .lt. 1 ) then ainc = -96. ak = -1. t63 = -1. a95 = -1. write (0,*) " ERROR: Inclination data missing in ARALEV" return end if if ( n .eq. 1 ) then ainc = xinc(1) ak = -1. t63 = t63max a95 = a95max ierr = 0 write (0,*) " WARNING: Only one observed inclination in ARALEV" return end if if ( n .gt. 10000 ) then ainc = -97. ak = -1. t63 = -1. a95 = -1. write (0,*) " ERROR: Too small dimension in ARALEV" return end if c Check if incl are out of range do i=1, n if ( xinc(i) .ge. -90.D0 .and. xinc(i) .le. 90.D0 ) goto 19 ainc = -98. ak = -1. t63 = -1. a95 = -1. write (0,*) " ERROR: Incl out of range [-90, +90] in ARALEV" return 19 continue end do c Check if all incl are identical do i=2, n if ( xinc(i) .ne. xinc(1) ) goto 29 end do ainc = xinc(1) ak = 1.E10 t63 = 0. a95 = 0. ierr = 0 write (0,*) " WARNING: All incl identical in ARALEV" return c Inclinations to co-inclinations 29 do i=1, n th(i) = 90.D0 - dble(xinc(i)) enddo c Calculate arithmetic mean to use as first guess s = 0.D0 s2 = 0.D0 c = 0.D0 do i=1, n s = s + th(i) s2 = s2 + th(i)**2 c = c + dcos(th(i)*dr) enddo c = c/dn rt = s/dn x = (s2 - s*s/dn)*dr*dr rk = 1.D10 if ( x/(dn-1.D0) .gt. 1.D-10 ) rk = (dn - 1.D0)/x rt1 = rt rk1 = rk c Iterate in the interior to find solution to (theta, kappa) c Start iteration at arithmetic mean (theta, kappa) rt = rt1 rk = rk1 ie1 = 0 the1 = rt akap1 = rk do j=1, 10000 rt = AL1 ( th, n, rt, rk ) rk = AL2 ( th, n, rt, rk ) dt = dabs(rt-the1) dk = dabs((rk-akap1)/rk) if ( j. gt. 10 .and. dt .lt. 1.D-6 1 .and. dk .lt. 1.D-6 ) goto 39 the1 = rt akap1 = rk enddo ie1 = 1 39 the1 = rt akap1 = rk xl1 = XLIK ( th, n, rt, rk ) c Find the maximum on the edge (theta = 0) rt = 0.D0 rk = rk1 akap2 = rk ie2 = 0 do j=1, 10000 x = COTH(rk) - c if ( x .gt. 1.D-10 ) then rk = 1.D0/x else rk = 1.D10 endif dk = dabs((rk-akap2)/rk) if ( j. gt. 4 .and. dk .lt. 1.D-6 ) goto 49 if ( rk .lt. 1.D-6 ) goto 49 akap2 = rk enddo ie2 = 1 49 the2 = 0.D0 akap2 = rk xl2 = XLIK ( th, n, rt, rk ) c Find the maximum on the edge (theta = 180) rt = 180.D0 rk = rk1 akap3 = rk ie3 = 0 do j=1, 10000 x = COTH(rk) + c if ( x .gt. 1.D-10 ) then rk = 1.D0/x else rk = 1.D10 endif dk = dabs((rk-akap3)/rk) if ( j. gt. 4 .and. dk .lt. 1.D-6 ) goto 59 if ( rk .lt. 1.D-6 ) goto 59 akap3 = rk enddo ie3 = 1 59 the3 = 180.D0 akap3 = rk xl3 = XLIK ( th, n, rt, rk ) c Find the maximum on the edge (kappa = 0) rt = 90.D0 rk = 0.D0 the4 = rt akap4 = rk xl4 = XLIK ( th, n, rt, rk ) c Use the best solution of the four isol = 1 ierr = ie1 if ( xl2 .gt. xl1 ) then the1 = the2 akap1 = akap2 xl1 = xl2 isol = 2 ierr = 1 if ( ie2 .eq. 0 ) ierr = 0 endif if ( xl3 .gt. xl1 ) then the1 = the3 akap1 = akap3 xl1 = xl3 isol = 3 ierr = 1 if ( ie3 .eq. 0 ) ierr = 0 endif if ( xl4 .gt. xl1 ) then the1 = the4 akap1 = akap4 xl1 = xl4 isol = 4 ierr = 0 endif ainc = 90.D0 - the1 ak = akap1 if ( ierr .ne. 0 ) 1 write (0,*) " WARNING: Convergence problems in ARALEV" c Test robustness of solution theta +/- 0.01° and kappa +/- 0.1% do i=1, 16 x = i rt = the1 + 0.01D0 * dcos( 22.5D0*x*dr ) if ( rt .lt. 0.D0 .or. rt .gt. 180D0 ) goto 69 rk = akap1 * ( 1.D0 + 0.001D0 * dsin( 22.5D0*x*dr ) ) xl = XLIK ( th, n, rt, rk ) if ( xl .gt. xl1 ) then ierr = ierr + 2 write (0,*) " WARNING: Robustness failure in ARALEV" endif 69 continue enddo c Estimation of angular standard deviation c Theta-63 calculated from (kappa), same method as Kono (1980) if ( akap1 .ge. 20.D0 ) co = 1.D0 + dlog(1.D0-0.63D0)/akap1 if ( akap1 .gt. 0.1D0 .and. akap1 .lt. 20.D0 ) 1 co = 1.D0 + dlog(1.D0 - 0.63D0* 2 (1.D0-dexp(-2.D0*akap1)))/akap1 if ( akap1 .le. 0.1D0 ) co = -0.26D0 + 0.4662D0*akap1 t63 = 90.D0 - dsign(90.D0,co) if ( dabs(co) .lt. 1.D0 ) 1 t63 = 90.D0 - datan( co/dsqrt(1.D0-co**2) )/dr if ( dble(t63) .gt. t63max ) t63 = t63max c Estimation of 95% (circular) symmetric confidence limit of the mean c Alpha-95 calculated from (N, kappa), same method as Kono (1980) co = 1.D0 - (dn - 1.D0) * (20.D0**(1.D0/(dn-1.D0)) - 1.D0) 1 /(dn*(akap1-1.D0) + 1.D0) a95 = 90.D0 - dsign(90.D0,co) if ( dabs(co) .lt. 1.D0 ) 1 a95 = 90.D0 - datan( co/dsqrt(1.D0-co**2) )/dr if ( dble(a95) .gt. a95max ) a95 = a95max return end c############################################################################ c c A L 1 - The Arason-Levi MLE Iteration Formula 1. c The method is described in a report and at: c Arason, P. and S. Levi (2006), The Maximum Likelihood Solution for c Inclination-only Data, Eos Trans. AGU, 87(52), Fall Meet. Suppl., c Abstract GP21B-1312. c See also: http://www.vedur.is/~arason/paleomag c c th input Vector of co-inclination data (in degrees) c n input Number of data in the sample c the input Current estimate of co-inclination (Theta) (in degrees) c ak input Current estimate of precision parameter (Kappa) c AL1 output New (better) estimate of Theta c c Author: c Dr. Pordur Arason, arason@vedur.is c Icelandic Meteorological Office, Bustadavegi 9, c IS-150 Reykjavik, Iceland c c History: c July 2006 Implementation and testing c Nov 2006 Comment clarifications c c############################################################################ real*8 function AL1 ( th, n, the, ak ) real*8 th(10000), the, ak integer n real*8 dr, s, c, x, bi0e, bi1e, bi1i0 dr = 0.0174532925199433D0 ! Degrees to radians (pi/180) s = 0.D0 c = 0.D0 do i=1, n x = ak*dsin(the*dr)*dsin(th(i)*dr) call BESSEL ( x, bi0e, bi1e, bi1i0 ) s = s + dsin(th(i)*dr) * bi1i0 c = c + dcos(th(i)*dr) enddo AL1 = datan2( s, c )/dr if ( AL1 .lt. 0.000001D0 ) AL1 = 0.000001D0 if ( AL1 .gt. 179.999999D0 ) AL1 = 179.999999D0 return end c############################################################################ c c A L 2 - The Arason-Levi MLE Iteration Formula 2. c The method is described in a report and at: c Arason, P. and S. Levi (2006), The Maximum Likelihood Solution for c Inclination-only Data, Eos Trans. AGU, 87(52), Fall Meet. Suppl., c Abstract GP21B-1312. c See also: http://www.vedur.is/~arason/paleomag c c th input Vector of co-inclination data (in degrees) c n input Number of data in the sample c the input Current estimate of co-inclination (Theta) (in degrees) c ak input Current estimate of precision parameter (Kappa) c AL2 output New (better) estimate of Kappa c c Author: c Dr. Pordur Arason, arason@vedur.is c Icelandic Meteorological Office, Bustadavegi 9, c IS-150 Reykjavik, Iceland c c History: c July 2006 Implementation and testing c Nov 2006 Comment clarifications c c############################################################################ real*8 function AL2 ( th, n, the, ak ) real*8 th(10000), the, ak integer n real*8 dr, s, c, x, bi0e, bi1e, bi1i0 real*8 COTH dr = 0.0174532925199433D0 ! Degrees to radians (pi/180) dn = n s = 0.D0 c = 0.D0 do i=1, n x = ak*dsin(the*dr)*dsin(th(i)*dr) call BESSEL ( x, bi0e, bi1e, bi1i0 ) s = s + dsin(th(i)*dr) * bi1i0 c = c + dcos(th(i)*dr) enddo x = dn*COTH(ak) - dcos(the*dr)*c - dsin(the*dr)*s AL2 = 1.D10 if ( x/dn .gt. 1.E-10) AL2 = dn/x if ( AL2 .lt. 1.D-06 ) AL2 = 1.D-06 return end c############################################################################ c c B E S S E L - Evaluation of the Hyperbolic Bessel functions c I0(x), I1(x) and their ratio I1(x)/I0(x). These functions are c sometimes also called the modified Bessel functions of order zero c and one. Since both the functions I0(x) and I1(x) increase c exponentially as x increases and become numerically unstable, c the output is given as I0(x)/exp(|x|), etc. c The ratio I1(x)/I0(x) increases smoothly from 0 to 1 as x c increases from zero. c c x input Range: -Inf < x < +Inf c bi0e output I0(x)/exp(|x|) Range: 0 < bi0e <= 1 c bi1e output I1(x)/exp(|x|) Range: -0.22 < bi1e < 0.22 c bi1i0 output I1(x)/I0(x) Range: -1 < bi1i0 < 1 c c Using the approximations of: c Press et al. (1989), Numerical recipes, The art of scientific c computing (Fortran version), Cambridge Univ. Press, Cambridge, 702 pp. c Abramowitz and Stegun (1972, eq. 9.8.1-4), Handbook of Mathematical c Functions with Formulas, Graphs, and Mathematical Tables, Dover, c New York. c c The accuracy of these approximations are: c abs[ exp(x)*bi0e - I0(x) ] < 1.6E-7 (for x<3.75) c abs[ exp(x)*bi1e - I1(x) ] < 0.3E-7 (for x<3.75) c abs[ sqrt(x)*(bi0e - I0(x)/exp(x)) ] < 1.9E-7 (for x>=3.75) c abs[ sqrt(x)*(bi1e - I1(x)/exp(x)) ] < 2.2E-7 (for x>=3.75) c c Author: c Dr. Pordur Arason, arason@vedur.is c Icelandic Meteorological Office, Reykjavik, Iceland c c History: c July 2006 Initial implementation c Nov 2006 Comment clarifications c c############################################################################ subroutine BESSEL ( x, bi0e, bi1e, bi1i0 ) real*8 x, bi0e, bi1e, bi1i0 real*8 t, b0 ,b1 real*8 p1, p2, p3, p4, p5, p6, p7 real*8 q1, q2, q3, q4, q5, q6, q7, q8, q9 real*8 u1, u2, u3, u4, u5, u6, u7 real*8 v1, v2, v3, v4, v5, v6, v7, v8, v9 data p1, p2, p3, p4, p5, p6, p7 / 1 1.D0, 3.5156229D0, 3.0899424D0, 1.2067492D0, 2 0.2659732D0, 0.360768D-1, 0.45813D-2 / data q1, q2, q3, q4, q5, q6, q7 ,q8, q9 / 1 0.39894228D0, 0.1328592D-1, 0.225319D-2, 2 -0.157565D-2, 0.916281D-2, -0.2057706D-1, 3 0.2635537D-1, -0.1647633D-1, 0.392377D-2 / data u1, u2, u3, u4, u5, u6, u7 / 1 0.5D0, 0.87890594D0, 0.51498869D0, 0.15084934D0, 2 0.2658733D-1, 0.301532D-2, 0.32411D-3 / data v1, v2, v3, v4, v5, v6, v7, v8, v9 / 1 0.39894228D0, -0.3988024D-1, -0.362018D-2, 2 0.163801D-2, -0.1031555D-1, 0.2282967D-1, 3 -0.2895312D-1, 0.1787654D-1, -0.420059D-2 / if ( dabs(x) .lt. 3.75D0 ) then t = (x/3.75D0)**2 b0 = p1+t*(p2+t*(p3+t*(p4+t*(p5+t*(p6+t*p7))))) b1 = x*(u1+t*(u2+t*(u3+t*(u4+t*(u5+t*(u6+t*u7)))))) bi0e = b0/dexp(dabs(x)) bi1e = b1/dexp(dabs(x)) bi1i0 = b1/b0 else t = 3.75D0/dabs(x) b0 = q1+t*(q2+t*(q3+t*(q4+t*(q5+t*(q6+t*(q7 1 + t*(q8+t*q9))))))) b1 = v1+t*(v2+t*(v3+t*(v4+t*(v5+t*(v6+t*(v7 1 + t*(v8+t*v9))))))) if ( x .lt. 0.D0 ) b1 = -b1 bi0e = b0/dsqrt(dabs(x)) bi1e = b1/dsqrt(dabs(x)) bi1i0 = b1/b0 endif return end c############################################################################ c c C O T H - Evaluation of the Hyperbolic Cotangens function coth(x). c c x input Range: -Inf < x < +Inf c x can not be 0 c c The function is similar to 1/x close to zero, and practically c identical to 1 for x>3. For the value x=0 the function returns c the value zero, although it should return +/-Inf. c We have the relation coth(-x) = - coth(x) c c The accuracy of these calculations are: c abs[ error ] < 1.E-13 (for 0.0001