c begin ellipse.for c------ ellipse.for j. c. lahr 4/6/84 c------ based on gpp3 computations c------ this program will plot one error ellipse c c dimension iaz(2),idp(2),se(3) c call pltts(0.,0.,0) call pltt(1.,1.,-3) c c input c c xzero x position of center of ellipse c yzero y position of center of ellipse c sinpk inches per km c iaz azimuth of first two ellipsoid axes c idp dip of first two ellipsoid axes c (the third axis is computed by cross multiplication) c se standard error of all three ellipsoid axes c ratio a scale factor for ellipses. set equal to 1.515/1.87 = .81 c for the joint 2-dimensional 68% confidence region c set equal to 1 to simply project the outline of the entire ellipsoid c phi angle local meridian makes with y plot axis c rot rotate the ellipse clockwise by this amount c ppaz azimuth of the cross section plane c morse 0 for map projection c 1 for cross section projection c idbug .ge. 2 for trace of computations in for010.dat c .ge. 3 for listing in for011.dat of x,y values xzero = 3.0 yzero = 3.0 sinpk = 1. iaz(1) = 45 iaz(2) = 135 idp(1) = 0 idp(2) = 0 se(1) = 4. se(2) = 2. se(3) = 1. elfac = 1. phi = 0. ppaz = -45. morse = 0 idbug = 2 rot = 0. 11 continue xzero = raskk('xzero', xzero) yzero = raskk('yzero', yzero) sinpk = raskk('sinpk', sinpk) iaz(1) = iaskk('iaz(1)', iaz(1)) idp(1) = iaskk('idp(1)', idp(1)) se(1) = raskk('se(1)', se(1)) iaz(2) = iaskk('iaz(2)', iaz(2)) idp(2) = iaskk('idp(2)', idp(2)) se(2) = raskk('se(2)', se(2)) se(3) = raskk('se(3)', se(3)) elfac = raskk('elfac', elfac) phi = raskk('phi', phi) morse = iaskk('map (0) or section (1)', morse) idbug = iaskk('idbug', idbug) rot = raskk('rot', rot) c icmout = 0 rpd = 1.745329251994e-2 c c c------ convert phi and ppaz to radians phi = phi*rpd ppaz = ppaz*rpd c call plttel(xzero,yzero,sinpk, 1 iaz,idp,se,ratio,phi,rot,ppaz,morse,idbug) go to 11 c 999 stop end subroutine plotel(xzero,yzero,sinpk, 1 iaz,idp,se,ratio,phi,rot,ppaz,morse,idbug) c c------- plot one standard deviation error ellipse. ------------------- c c c input c c xzero x position of center of ellipse c yzero y position of center of ellipse c sinpk inches per km c iaz azimuth of first two ellipsoid axes c idp dip of first two ellipsoid axes c (the third axis is computed by cross multiplication) c se standard error of all three ellipsoid axes c ratio a scale factor for ellipses. set equal to 1.515/1.87 = .81 c for the joint 2-dimensional 68% confidence region c set equal to 1 to simply project the outline of the entire ellipsoid c phi angle local meridian makes with y plot axis c rot rotate the ellipse clockwise by this amount c ppaz azimuth of the cross section plain c morse 0 for map projection c 1 for cross section projection c idbug .ge. 2 for trace of computations c output c dimension iaz(2),idp(2),se(3) c rpd = 1.745329251994e-2 c call ellips(iaz,idp,se,ratio,phi,ppaz,morse,idbug, 1 ael,bel,cel) c c------- the ellipse equation is: c 1 = ael*x**2 + bel*x*y + cel*y**2 c or in terms of r and theta: c 1 = ael*r**2*cos(theta)**2 + bel*r**2*cos(theta)*sin(theta) + c cel*r**2*sin(theta)**2 c c c------- plot a grid amax = 1./sqrt(ael) if(1./sqrt(cel) .gt. amax) amax = 1./sqrt(cel) call pltt(0.,0.,-998) imax = amax + 1.5 xi = -imax + xzero call pltt(xi, yzero, 3) do 100 i = -imax, imax xi = i + xzero call pltt(xi, yzero, 2) call pltt(xi, yzero+.1, 2) call pltt(xi, yzero, 2) 100 continue yi = -imax + yzero call pltt(xzero, yi, 3) do 200 i = -imax, imax yi = i + yzero call pltt(xzero, yi, 2) call pltt(xzero+.1, yi, 2) call pltt(xzero, yi, 2) 200 continue if(ael .eq. cel) go to 650 theta = 0.5*atan(bel/(ael-cel)) go to 675 650 theta = 45.0*rpd 675 dthet = 18.0*rpd do 700 i = 1,21 cs = cos(theta) sn = sin(theta) ra = 1.0/sqrt(ael*cs**2+bel*cs*sn+cel*sn**2) x = ra*cos(theta-rpd*rot)*sinpk + xzero y = ra*sin(theta-rpd*rot)*sinpk + yzero if(i .eq. 1) then if(idbug .ge. 3) write(11, *) 'p0 -999 0 0' call pltt(x,y,3) else if(idbug .ge. 3) write(11, *) 'p0 ', x, y, ' 0' call pltt(x,y,2) endif 700 theta = theta + dthet call pltt(x,y,3) call pltt(0.,0.,-998) return end subroutine ellips(iaz,idp,se,ratio,phi,ppaz,morse,idbug, 1 ael,bel,cel) c------- program to find shadow projection of the error ellipse on c a horizontal plane or on any vertical plane. c c input c c iaz azimuth of first two ellipsoid axes (degrees) c idp dip of first two ellipsoid axes (degrees) c se standard error of all three ellipsoid axes (kilometers) c ratio a scale factor for ellipses. set equal to 1.515/1.87 = .81 c for the joint 2-dimensional 68% confidence region c set equal to 1 to simply project the outline of the entire ellipsoid c phi angle the local meridian makes with y plot axis (radians) c ppaz azimuth of the cross section plain (radians) c morse 0 for map projection c 1 for cross section projection c idbug .ge. 2 for trace of computations c output c c the resulting ellipse is> ael*x**2 + bel*x*y + cel*y**2 = 1 c c c dimension iaz(2),idp(2),se(3) dimension pauv(3,3),tpacs(3,3),teacs(3,3) c data rpd/.017453292/ nofile = 10 c if(idbug .ge. 2) write(nofile,900) 900 format(/,' *** subroutine ellipse *** ') 50 do 150 i = 1,3 do 150 j = 1,3 tpacs(i,j) = 0.0 pauv(i,j) = 0.0 150 teacs(i,j) = 0.0 c c------- find ellipse principal axis vectors in terms of the projection c plane coordinates. z of the projection plane coords is the c projection direction. if(morse .eq. 1) go to 230 c c------- for map projection. ------------------------------------------- do 225 i = 1,2 dp = idp(i)*rpd az = iaz(i)*rpd app = az + phi if(idbug .le. 1) go to 220 write(nofile,1000) iaz(i),idp(i) 1000 format(' az = ',i3,7x,' dp = ',i3,7x) 220 pauv(i,1) = cos(dp)*sin(app) pauv(i,2) = cos(dp)*cos(app) 225 pauv(i,3) = -sin(dp) go to 240 c c------- for x-section projection. ------------------------------------- 230 do 235 i = 1,2 dp = idp(i)*rpd az = iaz(i)*rpd app = az + phi if(idbug .le. 1) go to 233 write(nofile,1000) iaz(i),idp(i) 233 pauv(i,1) = cos(dp)*cos(app-ppaz) pauv(i,2) = -sin(dp) 235 pauv(i,3) = cos(dp)*sin(app-ppaz) c c------- cross multiply to find third unit vector. 240 pauv(3,1) = pauv(1,2)*pauv(2,3)-pauv(2,2)*pauv(1,3) pauv(3,2) = pauv(2,1)*pauv(1,3)-pauv(1,1)*pauv(2,3) pauv(3,3) = pauv(1,1)*pauv(2,2)-pauv(2,1)*pauv(1,2) c if(idbug .le. 1) go to 250 write(nofile,1015) 1015 format(/1x,'transformation tensor from principal axis to ', 1 'projection plane coordinates'/) do 260 i = 1,3 260 write(nofile,1020) (pauv(i,j),j=1,3) 1020 format(3f12.2) c c------- set up tensor representing ellipsoid in principal coordinates. 250 do 275 i = 1,3 275 tpacs(i,i) = 1.0/(.0000001+(se(i))**2) c c------- transform this tensor to the projection coordinate system. do 375 i = 1,3 do 375 j = i,3 do 350 k = 1,3 do 350 l = 1,3 teacs(i,j) = teacs(i,j) + pauv(k,i)*pauv(l,j)*tpacs(k,l) 350 continue 375 teacs(j,i) = teacs(i,j) if(idbug .le. 1) go to 400 write(nofile,1050) 1050 format(/' error ellipsoid in principal coordinates', 1 ' and projection coordinates'/) do 380 i = 1,3 380 write(nofile,1075) (tpacs(i,j),j=1,3),(teacs(i,j),j=1,3) 1075 format(3f15.2,10x,3f15.2) c c------- ellipsoid is x.teacs.x - 1 = f = 0 c------- grad f is a normal vector. c------- grad f .k = 0 defines the plane> c------- teacs(1,3)*x + teacs(2,3)*y + teacs(3,3)*z = 0 c------- solve for z and substitute in f equation to obtain ellipse. --- c------- the ellipse is> ael*x**2 + bel*x*y + cel*y**2 = 1 c c------- sqrt(chi-square) for 3 degrees of freedom and 68% is 1.87 c sqrt(chi-square) for 2 degrees of freedom and 68% is 1.515 c ratio = 1.515/1.87 = .81 c ratio = 1.515/1.87 c a ratio of 0.81 scales the ellipse for two degrees of freedom, maintaining 68% c a ratio of 1.0 just projects the entire ellipsoid c ratio = 1. c 400 ael = ratio*(teacs(1,1)-teacs(1,3)*teacs(1,3)/teacs(3,3)) bel = ratio*2.0*(teacs(1,2)-teacs(1,3)*teacs(2,3)/teacs(3,3)) cel = ratio*(teacs(2,2)-teacs(2,3)*teacs(2,3)/teacs(3,3)) if(idbug .le. 1) return write(nofile,1100) ael,bel,cel 1100 format(/' ael =',e10.4,' bel =',e10.4,' cel =',e10.4) return end c end ellipse