program lavamax1 c performs maximum likelihood fit of 3D gaussian distribution c to paleomagnetic data read from external file parameter(nmax=20000,nmax2=20000) dimension rf(nmax2),ri(nmax2),rd(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(4) character*8 rlabel8 character*37 dumb character*20 rmodel character*80 string dimension pp(5,4),yy(5),ppp(4,4) external funcf,funci,funcfi,funcid,funcfid external gasdev common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid data pi/3.14159265358979323846264338328/ fac = pi/180. write(*,*)' ' write(*,*)'This program fits paleomagnetic vector data with a' write(*,*)'maximum-likelihood procedure, as outlined in' write(*,*)'Love, J. J. & Constable, C. G., 2003.' write(*,*)'Gaussian statistics for palaeomagnetic vectors,' write(*,*)'Geophys. J. Int., 152, 515-565.' write(*,*)' ' write(*,*)'The bulk of this program was written by J. J. Love.' write(*,*)'Several Numerical Recipes routines are used for' write(*,*)'performing the maximization of the likelihood.' write(*,*)' ' write(*,*)'The data format is like that generated by a companion' write(*,*)'program called lavadata.f.' write(*,*)' ' write(*,*)'Some of the subroutines used here are based on the' write(*,*)'partial summation of infinite series. If more' write(*,*)'precisionis needed, you can change nmax in rfdisti' write(*,*)'and f11 from 20 or so to a larger number ... ' write(*,*)'say 30 or so, but this, then, will probably require' write(*,*)'compilation in double precision.' write(*,*)' ' write(*,*)'The data must be entered in list form with fid' write(*,*)'(intensity, inclination, declination) triplets' write(*,*)'corresponding to coincident data. So, for example,' write(*,*)'separate intensity measurements and directional' write(*,*)'measurements which you wish to be associated,' write(*,*)'because they come (say) from the same sample or' write(*,*)'flow, should be entered on the same line. This,' write(*,*)'then, corresponds to a full vector triplet. On' write(*,*)'the other hand, if no intensity is available,' write(*,*)'but directions are, then only an inclination and' write(*,*)'declination are entered on a given line, with the' write(*,*)'missing intensity entered as a 999.9. Correspondignly,' write(*,*)'missing inclinations are entered as 99.9, and' write(*,*)'missing declinations are entered as 999.9. The' write(*,*)'program will count all of the different data groupings' write(*,*)'and then automatically use the appropriate pdf from' write(*,*)'Love & Constable 2003 to execute a maximum-likelihood' write(*,*)'estimate.' write(*,*)' ' write(*,*)'Enter the data file name now.' write(*,*)' ' read(*,'(a20)')rmodel write(string,'(a20)')rmodel open(1,file=string,status='unknown') 8989 format(a20) write(*,*)' ' write(*,*)'Data will be normalized to common latitude.' write(*,*)' ' write(*,*)'Enter that latitude now.' write(*,*)' ' read(*,*)rlatfix rincfix = atan(2*tan(rlatfix*fac)) ffactorfix = sqrt(1+3.*sin(rlatfix*fac)**2) write(*,*)' ' write(*,*)'Now that you have entered the data you need to exploit' write(*,*)'them. What kind of fit do you wish to perform?' write(*,*)' ' write(*,*)'An f-fit is one where only intensities are fitted,' write(*,*)'even though there may be some associated directions,' write(*,*)'they will not be used in the fitting procedure.' write(*,*)' ' write(*,*)'An i-fit is one where only inclinations are fitted.' write(*,*)' ' write(*,*)'etc.' write(*,*)' ' write(*,*)'An fi-fit is one where intensity-inclination' write(*,*)'pairs, plus any separate intensity and inclination' write(*,*)'data, are fitted. ' write(*,*)' ' write(*,*)'An id-fit is one where inclination-declination pairs,' write(*,*)'plus any separate inclinations that may come from' write(*,*)'azimuthally unoriented samples, are fitted. ' write(*,*)' ' write(*,*)'An fid-fit is one where intensity-inclination-' write(*,*)'declination triplets, plus any intensity-inclination' write(*,*)'pairs, any inclination-declination pairs, and any' write(*,*)'separate intensities and inclinations, are fitted.' write(*,*)'This last option, of course, makes maximum' write(*,*)'exploitation of the informational content of' write(*,*)'the data.' write(*,*)' ' write(*,*)'Enter 1 f fit' write(*,*)'Enter 2 i fit' write(*,*)'Enter 3 fi fit' write(*,*)'Enter 4 id fit' write(*,*)'Enter 5 fid fit' write(*,*)' ' read(*,*)ifit write(*,*)' ' write(*,*)'Now you need to set the convergence tolerance.' write(*,*)'This should be a small number like 1e-7 or so.' write(*,*)'Smaller numbers will yield more accurate results,' write(*,*)'but will take longer to obtain.' write(*,*)' ' write(*,*)'Enter the tolerance now.' write(*,*)' ' read(*,*)ftol c load in the data do i = 1,nmax2 read(1,'(a37,f5.1,f6.1,a23)',end=2222)dumb,theta,phi rinc1 = atan(2*tan(theta*fac)) ffactor1 = sqrt(1+3.*sin(theta*fac)**2) do j = 1,nmax2 read(1,'(a8,i3,f6.1,1x,f5.1,1x,f6.1,1x,f6.1,1x,f6.1)') + rlabel8,nnn,rf1,rfsig1,rd1,ri1,alpha if (rlabel8 .eq. '1000 ') then go to 1111 end if if (rd1 .gt. 180. .and. rd1 .ne. 999.9) rd1 = rd1 - 360 c fix data for latitude if (rf1 .ne. 999.9 .and. ri1 .eq. 99.9 .and. + rd1 .eq. 999.9) then rf1 = rf1*ffactorfix/ffactor1 nf = nf + 1 rf(nf) = rf1 rff = rff + rf1 mff = mff + 1 else if (rf1 .eq. 999.9 .and. ri1 .ne. 99.9 .and. + rd1 .eq. 999.9) then ri1 = ri1 + (rincfix - rinc1)/fac ni = ni + 1 ri(ni) = ri1*fac else if (rf1 .ne. 999.9 .and. ri1 .ne. 99.9 .and. + rd1 .eq. 999.9) then rf1 = rf1*ffactorfix/ffactor1 ri1 = ri1 + (rincfix - rinc1)/fac nfi = nfi + 1 rfi(1,nfi) = rf1 rfi(2,nfi) = ri1*fac rff = rff + rf1 mff = mff + 1 else if (rf1 .eq. 999.9 .and. ri1 .ne. 99.9 .and. + rd1 .ne. 999.9) then ri1 = ri1 + (rincfix - rinc1)/fac nid = nid + 1 rid(1,nid) = ri1*fac rid(2,nid) = rd1*fac rx1 = rx1 + cos(ri1*fac)*cos(rd1*fac) ry1 = ry1 + cos(ri1*fac)*sin(rd1*fac) rz1 = rz1 + sin(ri1*fac) mdd = mdd + 1 else if (rf1 .ne. 999.9 .and. ri1 .ne. 99.9 .and. + rd1 .ne. 999.9) then rf1 = rf1*ffactorfix/ffactor1 ri1 = ri1 + (rincfix - rinc1)/fac nfid = nfid + 1 rfid(1,nfid) = rf1 rfid(2,nfid) = ri1*fac rfid(3,nfid) = rd1*fac rx1 = rx1 + cos(ri1*fac)*cos(rd1*fac) ry1 = ry1 + cos(ri1*fac)*sin(rd1*fac) rz1 = rz1 + sin(ri1*fac) rff = rff + rf1 mdd = mdd + 1 mff = mff + 1 end if end do 1111 continue end do 2222 continue c do the appropriate fit if (ifit .eq. 1) then write(*,*)' ' write(*,*)'Enter a guess for f (microT)' write(*,*)'and the dispersion (sigma) (microT).' write(*,*)' ' read(*,*)rmfsave,rmssave write(*,*)' ' write(*,*)'Computer is working now ...' dumb1 = rmfsave dumb2 = rmssave pp(1,1) = dumb1 + gasdev(idumb)*0.01 pp(1,2) = dumb2 + gasdev(idumb)*0.01 pp(2,1) = dumb1 + gasdev(idumb)*0.01 pp(2,2) = dumb2 + gasdev(idumb)*0.01 pp(3,1) = dumb1 + gasdev(idumb)*0.01 pp(3,2) = dumb2 + gasdev(idumb)*0.01 do i = 1,3 do j = 1,2 p(j) = pp(i,j) end do yy(i) = funcf(p) end do call amoeba(pp,yy,5,4,2,ftol,funcf,iter) x = sqrt((pp(1,1)**2)/3.) y = x z = x s = pp(1,2) write(*,*)' ' write(*,*)'Maximum-likelihood results for f-fit:' write(*,*)' ' write(*,*)'Cartesian components are generic, since directional' write(*,*)'data were not used.' write(*,*)' ' write(*,*)' x y z sigma' write(*,*)' ' write(*,'(4(1x,f9.4))')x,y,z,s write(*,*)' ' fff = sqrt(x**2+y**2+z**2) ssf = s/fff write(*,*)' f sig/f' write(*,*)' ' write(*,'(6(1x,f9.4))')fff,ssf write(*,*)' ' else if (ifit .eq. 2) then write(*,*)' ' write(*,*)'Enter a guess for i (in degrees) and the' write(*,*)'relative dispersion (sig/f).' write(*,*)' ' read(*,*)rmisave,rmfssave rmisave = rmisave*fac write(*,*)' ' write(*,*)'Computer is working now ... ' write(*,*)'... and this may take some time.' dumb1 = 1./rmfssave dumb2 = rmisave pp(1,1) = dumb1 + gasdev(idumb)*0.01 pp(1,2) = dumb2 + gasdev(idumb)*0.01 pp(2,1) = dumb1 + gasdev(idumb)*0.01 pp(2,2) = dumb2 + gasdev(idumb)*0.01 pp(3,1) = dumb1 + gasdev(idumb)*0.01 pp(3,2) = dumb2 + gasdev(idumb)*0.01 do i = 1,3 do j = 1,2 p(j) = pp(i,j) end do yy(i) = funci(p) end do call amoeba(pp,yy,5,4,2,ftol,funci,iter) x = 30.*cos(pp(1,2)) y = 0 z = 30.*sin(pp(1,2)) s = 30./pp(1,1) write(*,*)' ' write(*,*)'Maximum-likelihood results for i-fit:' write(*,*)' ' write(*,*)'Cartesian components are generic, since' write(*,*)'intensity and declination data were not used.' write(*,*)' ' write(*,*)' x y z sigma' write(*,*)' ' write(*,'(4(1x,f9.4))')x,y,z,s write(*,*)' ' fff = sqrt(x**2+y**2+z**2) hhh = sqrt(x**2+y**2) rii = atan2(z,hhh)/fac ssf = s/fff write(*,*)' i sig/f' write(*,*)' ' write(*,'(6(1x,f9.4))')rii,ssf write(*,*)' ' else if (ifit .eq. 3) then write(*,*)' ' write(*,*)'Enter a guess for f (microT), i (degrees)' write(*,*)'and the dispersion (sigma) (microT).' write(*,*)' ' read(*,*)rmfsave,rmisave,rmssave write(*,*)' ' write(*,*)'Computer is working now ...' dumb1 = rmfsave dumb2 = rmisave*fac dumb3 = rmssave pp(1,1) = dumb1 + gasdev(idumb)*0.01 pp(1,2) = dumb2 + gasdev(idumb)*0.01 pp(1,3) = dumb3 + gasdev(idumb)*0.01 pp(2,1) = dumb1 + gasdev(idumb)*0.01 pp(2,2) = dumb2 + gasdev(idumb)*0.01 pp(2,3) = dumb3 + gasdev(idumb)*0.01 pp(3,1) = dumb1 + gasdev(idumb)*0.01 pp(3,2) = dumb2 + gasdev(idumb)*0.01 pp(3,3) = dumb3 + gasdev(idumb)*0.01 pp(4,1) = dumb1 + gasdev(idumb)*0.01 pp(4,2) = dumb2 + gasdev(idumb)*0.01 pp(4,3) = dumb3 + gasdev(idumb)*0.01 do i = 1,4 do j = 1,3 p(j) = pp(i,j) end do yy(i) = funcfi(p) end do call amoeba(pp,yy,5,4,3,ftol,funcfi,iter) x = pp(1,1)*cos(pp(1,2)) y = 0 z = pp(1,1)*sin(pp(1,2)) s = pp(1,3) write(*,*)' ' write(*,*)'Maximum-likelihood results:' write(*,*)' ' write(*,*)'Cartesian components are generic, since' write(*,*)'declination data were not used.' write(*,*)' ' write(*,*)' x y z sigma' write(*,*)' ' write(*,'(4(1x,f9.4))')x,y,z,s write(*,*)' ' fff = sqrt(x**2+y**2+z**2) hhh = sqrt(x**2+y**2) rii = atan2(z,hhh)/fac rdd = atan2(y,x)/fac ssf = s/fff write(*,*)' f h i sig/f' write(*,*)' ' write(*,'(6(1x,f9.4))')fff,hhh,rii,ssf write(*,*)' ' else if (ifit .eq. 4) then write(*,*)' ' write(*,*)'Enter a guess for i, d (in degrees) and' write(*,*)'the relative dispersion (sig/f).' write(*,*)' ' read(*,*)rmisave,rmdsave,rmfssave write(*,*)' ' write(*,*)'Computer is working now ...' dumb1 = 1./rmfssave dumb2 = rmisave*fac dumb3 = rmdsave*fac pp(1,1) = dumb1 + gasdev(idumb)*0.01 pp(1,2) = dumb2 + gasdev(idumb)*0.01 pp(1,3) = dumb3 + gasdev(idumb)*0.01 pp(2,1) = dumb1 + gasdev(idumb)*0.01 pp(2,2) = dumb2 + gasdev(idumb)*0.01 pp(2,3) = dumb3 + gasdev(idumb)*0.01 pp(3,1) = dumb1 + gasdev(idumb)*0.01 pp(3,2) = dumb2 + gasdev(idumb)*0.01 pp(3,3) = dumb3 + gasdev(idumb)*0.01 pp(4,1) = dumb1 + gasdev(idumb)*0.01 pp(4,2) = dumb2 + gasdev(idumb)*0.01 pp(4,3) = dumb3 + gasdev(idumb)*0.01 do i = 1,4 do j = 1,3 p(j) = pp(i,j) end do yy(i) = funcid(p) end do call amoeba(pp,yy,5,4,3,ftol,funcid,iter) x = 30.*cos(pp(1,2))*cos(pp(1,3)) y = 30.*cos(pp(1,2))*sin(pp(1,3)) z = 30.*sin(pp(1,2)) s = 30./pp(1,1) write(*,*)' ' write(*,*)'Maximum-likelihood results:' write(*,*)' ' write(*,*)'Cartesian components are generic, since' write(*,*)'intensity data were not used.' write(*,*)' ' write(*,*)' x y z sigma' write(*,*)' ' write(*,'(4(1x,f9.4))')x,y,z,s write(*,*)' ' fff = sqrt(x**2+y**2+z**2) hhh = sqrt(x**2+y**2) rii = atan2(z,hhh)/fac rdd = atan2(y,x)/fac ssf = s/fff write(*,*)' i d sig/f' write(*,*)' ' write(*,'(6(1x,f9.4))')rii,rdd,ssf write(*,*)' ' else if (ifit .eq. 5) then write(*,*)' ' write(*,*)'Enter a guess for f (microT) ,i,d (degrees)' write(*,*)'and the dispersion (sigma) (microT).' write(*,*)' ' read(*,*)rmfsave,rmisave,rmdsave,rmssave write(*,*)' ' write(*,*)'Computer is working now ...' dumb1 = rmfsave dumb2 = rmisave*fac dumb3 = rmdsave*fac dumb4 = rmssave pp(1,1) = dumb1 + gasdev(idumb)*0.01 pp(1,2) = dumb2 + gasdev(idumb)*0.01 pp(1,3) = dumb3 + gasdev(idumb)*0.01 pp(1,4) = dumb4 + gasdev(idumb)*0.01 pp(2,1) = dumb1 + gasdev(idumb)*0.01 pp(2,2) = dumb2 + gasdev(idumb)*0.01 pp(2,3) = dumb3 + gasdev(idumb)*0.01 pp(2,4) = dumb4 + gasdev(idumb)*0.01 pp(3,1) = dumb1 + gasdev(idumb)*0.01 pp(3,2) = dumb2 + gasdev(idumb)*0.01 pp(3,3) = dumb3 + gasdev(idumb)*0.01 pp(3,4) = dumb4 + gasdev(idumb)*0.01 pp(4,1) = dumb1 + gasdev(idumb)*0.01 pp(4,2) = dumb2 + gasdev(idumb)*0.01 pp(4,3) = dumb3 + gasdev(idumb)*0.01 pp(4,4) = dumb4 + gasdev(idumb)*0.01 pp(5,1) = dumb1 + gasdev(idumb)*0.01 pp(5,2) = dumb2 + gasdev(idumb)*0.01 pp(5,3) = dumb3 + gasdev(idumb)*0.01 pp(5,4) = dumb4 + gasdev(idumb)*0.01 do i = 1,5 do j = 1,4 p(j) = pp(i,j) end do yy(i) = funcfid(p) end do call amoeba(pp,yy,5,4,4,ftol,funcfid,iter) x = pp(1,1)*cos(pp(1,2))*cos(pp(1,3)) y = pp(1,1)*cos(pp(1,2))*sin(pp(1,3)) z = pp(1,1)*sin(pp(1,2)) s = pp(1,4) write(*,*)' ' write(*,*)'Maximum-likelihood results for fid-fit:' write(*,*)' ' write(*,*)' x y z sigma' write(*,*)' ' write(*,'(4(1x,f9.4))')x,y,z,s write(*,*)' ' fff = sqrt(x**2+y**2+z**2) hhh = sqrt(x**2+y**2) rii = atan2(z,hhh)/fac rdd = atan2(y,x)/fac ssf = s/fff write(*,*)' f h i d sig/f' write(*,*)' ' write(*,'(6(1x,f9.4))')fff,hhh,rii,rdd,ssf write(*,*)' ' end if stop end function funcf(p) parameter(nmax2=20000) dimension rf1(nmax2) dimension rf(nmax2),ri(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(*) common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid funcf = 0 rmf = p(1) rms = p(2) if (nf .ne. 0) then call ffunc(rmf,rms,nf,rf,rl) funcf = funcf - rl end if if (nfi .ne. 0) then do i = 1,nfi rf1(i) = rfi(1,i) end do call ffunc(rmf,rms,nfi,rf1,rl) funcf = funcf - rl end if if (nfid .ne. 0) then do i = 1,nfid rf1(i) = rfid(1,i) end do call ffunc(rmf,rms,nfid,rf1,rl) funcf = funcf - rl end if return end function funci(p) parameter(nmax2=20000) dimension ri1(nmax2) dimension rf(nmax2),ri(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(*) common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid data pi/3.14159265358979323846264338328/ funci = 0 rmfs = p(1) rmi = p(2) if (ni .ne. 0) then call ifunc(rmfs,rmi,ni,ri,rl) funci = funci - rl end if if (nfi .ne. 0) then do i = 1,nfi ri1(i) = rfi(2,i) end do call ifunc(rmfs,rmi,nfi,ri1,rl) funci = funci - rl end if if (nid .ne. 0) then do i = 1,nid ri1(i) = rid(1,i) end do call ifunc(rmfs,rmi,nid,ri1,rl) funci = funci - rl end if if (nfid .ne. 0) then do i = 1,nfid ri1(i) = rfid(2,i) end do call ifunc(rmfs,rmi,nfid,ri1,rl) funci = funci - rl end if return end function funcfi(p) parameter(nmax2=20000) dimension rfi1(2,nmax2),ri1(nmax2) dimension rf(nmax2),ri(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(*) common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid data pi/3.14159265358979323846264338328/ funcfi = 0 rmfs = p(1)/p(3) rmf = p(1) rmi = p(2) rms = p(3) if (nf .ne. 0) then call ffunc(rmf,rms,nf,rf,rl) funcfi = funcfi - rl end if if (ni .ne. 0) then call ifunc(rmfs,rmi,ni,ri,rl) funcfi = funcfi - rl end if if (nfi .ne. 0) then call fifunc(rmf,rmi,rms,nfi,rfi,rl) funcfi = funcfi - rl end if if (nid .ne. 0) then do i = 1,nid ri1(i) = rid(1,i) end do call ifunc(rmfs,rmi,nid,ri1,rl) funcfi = funcfi - rl end if if (nfid .ne. 0) then do i = 1,nfid rfi1(1,i) = rfid(1,i) rfi1(2,i) = rfid(2,i) end do call fifunc(rmf,rmi,rms,nfid,rfi1,rl) funcfi = funcfi - rl end if return end function funcid(p) parameter(nmax2=20000) dimension rid1(2,nmax2),ri1(nmax2) dimension rf(nmax2),ri(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(*) common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid data pi/3.14159265358979323846264338328/ funcid = 0 rmfs = p(1) rmi = p(2) rmd = p(3) if (ni .ne. 0) then call ifunc(rmfs,rmi,ni,ri,rl) funcid = funcid - rl end if if (nfi .ne. 0) then do i = 1,nfi ri1(i) = rfi(2,i) end do call ifunc(rmfs,rmi,nfi,ri1,rl) funcid = funcid - rl end if if (nid .ne. 0) then call idfunc(rmfs,rmi,rmd,nid,rid,rl) funcid = funcid - rl end if if (nfid .ne. 0) then do i = 1,nfid rid1(1,i) = rfid(2,i) rid1(2,i) = rfid(3,i) end do call idfunc(rmfs,rmi,rmd,nfid,rid1,rl) funcid = funcid - rl end if return end function funcfid(p) parameter(nmax2=20000) dimension rf(nmax2),ri(nmax2) dimension rfi(2,nmax2),rid(2,nmax2),rfid(3,nmax2) dimension p(*) common /data/ rf,ri,rfi,rid,rfid common /ndata/ nf,ni,nfi,nid,nfid data pi/3.14159265358979323846264338328/ funcfid = 0 rmfs = p(1)/p(4) rmf = p(1) rmi = p(2) rmd = p(3) rms = p(4) if (nf .ne. 0) then call ffunc(rmf,rms,nf,rf,rl) funcfid = funcfid - rl end if if (ni .ne. 0) then call ifunc(rmfs,rmi,ni,ri,rl) funcfid = funcfid - rl end if if (nfi .ne. 0) then call fifunc(rmf,rmi,rms,nfi,rfi,rl) funcfid = funcfid - rl end if if (nid .ne. 0) then call idfunc(rmfs,rmi,rmd,nid,rid,rl) funcfid = funcfid - rl end if if (nfid .ne. 0) then call fidfunc(rmf,rmi,rmd,rms,nfid,rfid,rl) funcfid = funcfid - rl end if return end subroutine ffunc(rmf,rms,n,rf,rl) dimension rf(*) rl = 0 do i = 1,n rf1 = rf(i) call rfdistf(rmf,rms,rf1,distf) rl = rl + log(distf) end do return end subroutine ifunc(rmfs,rmi,n,ri,rl) dimension ri(*) rl = 0 do i = 1,n ri1 = ri(i) call rfdisti(rmfs,rmi,ri1,disti) rl = rl + log(disti) end do return end subroutine fifunc(rmf,rmi,rms,n,rfi,rl) dimension rfi(2,*) rl = 0 do i = 1,n rf1 = rfi(1,i) ri1 = rfi(2,i) call rfdistfi(rmf,rmi,rms,rf1,ri1,distfi) rl = rl + log(distfi) end do return end subroutine idfunc(rmfs,rmi,rmd,n,rid,rl) dimension rid(2,*) rl = 0 do i = 1,n ri1 = rid(1,i) rd1 = rid(2,i) call rfdistid(rmfs,rmi,rmd,ri1,rd1,distid) rl = rl + log(distid) end do return end subroutine fidfunc(rmf,rmi,rmd,rms,n,rfid,rl) dimension rfid(3,*) rl = 0 do i = 1,n rf1 = rfid(1,i) ri1 = rfid(2,i) rd1 = rfid(3,i) call rfdistfid(rmf,rmi,rmd,rms,rf1,ri1,rd1, + distfid) rl = rl + log(distfid) end do return end subroutine rfdistfid(rmf,rmi,rmd,rms,rf,ri,rd, + distfid) data pi/3.14159265358979323846264338328/ c formula 11 distfid = 0 dumb = rf**2*cos(ri)/((2*pi)**(1.5)*rms**3)* + exp(-(rf*cos(ri)*cos(rd)-rmf*cos(rmi)*cos(rmd))**2/ + (2*rms**2))* + exp(-(rf*cos(ri)*sin(rd)-rmf*cos(rmi)*sin(rmd))**2/ + (2*rms**2))* + exp(-(rf*sin(ri)-rmf*sin(rmi))**2/ + (2*rms**2)) distfid = dumb return end subroutine rfdistf(rmf,rms,rf,distf) data pi/3.14159265358979323846264338328/ c formula D1 distf = sqrt(2./pi)*rf/(rms*rmf)*sinh(rf*rmf/rms**2)* + exp(-(rf**2)/(2*rms**2) - (rmf**2)/(2*rms**2)) return end subroutine rfdisti(rmfs,rmi,ri,disti) data pi/3.14159265358979323846264338328/ c formula F1 nmax = 20 c if more precision is needed, you can change nmax to a larger number c but this, then, will probably require compilation in double precision disti1 = 0 disti2 = 0 distic = 0 dumbc = rmfs*cos(rmi)*cos(ri) dumbs = rmfs*sin(rmi)*sin(ri) do m = 0,nmax m2 = 2*m do n = 0,nmax n2 = 2*n disti1 = disti1 + exp(factfactln(m2+n2+1) - + factfactln(m2) - factfactln(m2) - factln(n2))* + dumbc**m2*dumbs**n2 disti2 = disti2 + exp(factfactln(m2+n2+2) - + factfactln(m2) - factfactln(m2) - factln(n2+1))* + dumbc**m2*dumbs**(n2+1) end do end do dumb = cos(ri)*exp(-0.5*rmfs**2) disti = dumb*(disti1/2.+disti2/sqrt(2*pi)) return end subroutine rfdisti1(rmfs,rmi,ri,disti) data pi/3.14159265358979323846264338328/ c formula F1 nmax = 20 c if more precision is needed, you can change nmax to a larger number c but this, then, will probably require compilation in double precision disti1 = 0 disti2 = 0 dumb1 = 0.5*(rmfs*cos(rmi)*cos(ri))**2 dumb2 = rmfs*sin(rmi)*sin(ri) do n = 0,nmax n2 = 2*n disti1 = disti1 + (n2+1)/exp(factfactln(n2))* + dumb2**n2*f11(n+1.5,1.,dumb1) disti2 = disti2 + (n2+2)/exp(factfactln(n2+1))* + dumb2**(n2+1)*f11(n+2.,1.,dumb1) end do dumb = cos(ri)*exp(-0.5*rmfs**2) disti = dumb*(disti1/2.+disti2/sqrt(2*pi)) return end function f11(a,c,z) nmax = 20 c if more precision is needed, you can change nmax to a larger number c but this, then, will probably require compilation in double precision f11 = 1 do n = 0,nmax dumb = exp(gammln(a+n) - gammln(c+n) - factln(n))*z**n f11 = f11 + dumb end do f11 = f11*exp(gammln(c)-gammln(a)) return end subroutine rfdistfi(rmf,rmi,rms,rf,ri,distfi) data pi/3.14159265358979323846264338328/ c formula C1 dumbs = rmf*rf*sin(rmi)*sin(ri)/rms**2 dumbc = rmf*rf*cos(rmi)*cos(ri)/rms**2 distfi = bessi0(dumbc)*exp(dumbs) dumb = rf**2*cos(ri)* + exp(-0.5*(rf/rms)**2-0.5*(rmf/rms)**2)/(sqrt(2*pi)*rms**3) distfi = distfi*dumb return end FUNCTION bessi0(x) REAL bessi0,x REAL ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0, *1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, *0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1, *-0.1647633d-1,0.392377d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. function hypser(a,b,c,z) INTEGER n fac=1 temp=fac aa=a bb=b cc=c do 11 n=1,100 fac=fac*aa*bb/cc fac=fac*z/n hypser=temp+fac if (hypser.eq.temp) return temp=hypser aa=aa+1. bb=bb+1. cc=cc+1. 11 continue return pause 'convergence failure in hypser' END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. subroutine rfdistid(rmfs,rmi,rmd,ri,rd,distid) data pi/3.14159265358979323846264338328/ c formula A1 dumb = cos(rmi)*cos(rmd)*cos(ri)*cos(rd) + + cos(rmi)*sin(rmd)*cos(ri)*sin(rd) + sin(rmi)*sin(ri) distid = (1./(4*pi))*cos(ri)*exp(-0.5*rmfs**2)* + ((rmfs**2*dumb**2 + 1)*exp(0.5*rmfs**2*dumb**2)* + (1+erf(rmfs*dumb/sqrt(2.)))+sqrt(2./pi)*rms*dumb) return end FUNCTION factln(n) INTEGER n REAL factln CU USES gammln REAL a(100),gammln SAVE a DATA a/100*-1./ if (n.lt.0) then pause 'negative factorial in factln' end if if (n.le.99) then if (a(n+1).lt.0.) a(n+1)=gammln(n+1.) factln=a(n+1) else factln=gammln(n+1.) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION gammln(xx) REAL gammln,xx INTEGER j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, *-.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION gasdev(idum) INTEGER idum REAL gasdev CU USES ran2 INTEGER iset REAL fac,gset,rsq,v1,v2,ran2 SAVE iset,gset DATA iset/0/ if (iset.eq.0) then 1 v1=2.*ran2(idum)-1. v2=2.*ran2(idum)-1. rsq=v1**2+v2**2 if(rsq.ge.1..or.rsq.eq.0.)goto 1 fac=sqrt(-2.*log(rsq)/rsq) gset=v1*fac gasdev=v2*fac iset=1 else gasdev=gset iset=0 endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION ran2(idum) INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter) INTEGER iter,mp,ndim,np,NMAX,ITMAX REAL ftol,p(mp,np),y(mp),funk PARAMETER (NMAX=20,ITMAX=5000) EXTERNAL funk CU USES amotry,funk INTEGER i,ihi,ilo,inhi,j,m,n REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry c Simplex code taken from Numerical Recipes iter=0 1 do 12 n=1,ndim sum=0. do 11 m=1,ndim+1 sum=sum+p(m,n) 11 continue psum(n)=sum 12 continue 2 ilo=1 if (y(1).gt.y(2)) then ihi=1 inhi=2 else ihi=2 inhi=1 endif do 13 i=1,ndim+1 if(y(i).le.y(ilo)) ilo=i if(y(i).gt.y(ihi)) then inhi=ihi ihi=i else if(y(i).gt.y(inhi)) then if(i.ne.ihi) inhi=i endif 13 continue rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) if (rtol.lt.ftol) then swap=y(1) y(1)=y(ilo) y(ilo)=swap do 14 n=1,ndim swap=p(1,n) p(1,n)=p(ilo,n) p(ilo,n)=swap 14 continue return endif c if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba' iter=iter+2 ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0) if (ytry.le.y(ilo)) then ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0) else if (ytry.ge.y(inhi)) then ysave=y(ihi) ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5) if (ytry.ge.ysave) then do 16 i=1,ndim+1 if(i.ne.ilo)then do 15 j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j)) p(i,j)=psum(j) 15 continue y(i)=funk(psum) endif 16 continue iter=iter+ndim goto 1 endif else iter=iter-1 endif goto 2 END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac) INTEGER ihi,mp,ndim,np,NMAX REAL amotry,fac,p(mp,np),psum(np),y(mp),funk PARAMETER (NMAX=20) EXTERNAL funk CU USES funk INTEGER j REAL fac1,fac2,ytry,ptry(NMAX) fac1=(1.-fac)/ndim fac2=fac1-fac do 11 j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 11 continue ytry=funk(ptry) if (ytry.lt.y(ihi)) then y(ihi)=ytry do 12 j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) 12 continue endif amotry=ytry return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION erf(x) REAL erf,x CU USES gammp REAL gammp if(x.lt.0.)then erf=-gammp(.5,x**2) else erf=gammp(.5,x**2) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION gammp(a,x) REAL a,gammp,x CU USES gcf,gser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammp' if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammp=gamser else call gcf(gammcf,a,x,gln) gammp=1.-gammcf endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. SUBROUTINE gcf(gammcf,a,x,gln) INTEGER ITMAX REAL a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) CU USES gammln INTEGER i REAL an,b,c,d,del,h,gammln gln=gammln(a) b=x+1.-a c=1./FPMIN d=1./b h=d do 11 i=1,ITMAX an=-i*(i-a) b=b+2. d=an*d+b if(abs(d).lt.FPMIN)d=FPMIN c=b+an/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue pause 'a too large, ITMAX too small in gcf' 1 gammcf=exp(-x+a*log(x)-gln)*h return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. SUBROUTINE gser(gamser,a,x,gln) INTEGER ITMAX REAL a,gamser,gln,x,EPS PARAMETER (ITMAX=100,EPS=3.e-7) CU USES gammln INTEGER n REAL ap,del,sum,gammln gln=gammln(a) if(x.le.0.)then if(x.lt.0.)pause 'x < 0 in gser' gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,ITMAX ap=ap+1. del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*EPS)goto 1 11 continue pause 'a too large, ITMAX too small in gser' 1 gamser=sum*exp(-x+a*log(x)-gln) return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. FUNCTION factfactln(n) INTEGER n REAL factfactln CU USES gammln REAL a(100),gammln SAVE a data pi/3.14159265358979323846264338328/ DATA a/100*-1./ if (n.lt.-1) then pause 'negative factorial in factln' end if if (n.eq.-1) then factfactln=0 return end if if (mod(n,2) .eq. 0) then nn = n/2 if (n.le.99) then if (a(n+1).lt.0.) a(n+1)=factln(nn)+nn*log(2.) factfactln=a(n+1) else factfactln=factln(nn)+nn*log(2.) endif else rn = float(n)/2. if (n.le.99) then if (a(n+1).lt.0.) a(n+1)=gammln(rn+1.)+ + rn*log(2.)+0.5*log(2./pi) factfactln=a(n+1) else factfactln=gammln(rn+1.)+rn*log(2.)+0.5*log(2./pi) endif end if return END FUNCTION bessi1(x) REAL bessi1,x REAL ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0, *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1, *0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{.