program ArcXsect implicit none real, parameter :: PI = 3.141592653589793 real, parameter :: D2R = PI/180. real, parameter :: R2D = 180./PI type cartPnt real :: x, y, z end type cartPnt type geoditic real :: lon real :: lat end type geoditic type arc type(geoditic) :: sPnt type(cartPnt) :: sxyz type(geoditic) :: ePnt type(cartPnt) :: exyz end type arc integer :: n real :: mag, radius real :: tmp, tmpsum real :: lambda, theta real :: roundoff real :: Angles(3) real :: N1(3), N2(3), L(3), I1(3), XSECT(3) real :: Spnt(3), Epnt(3) logical :: found type(arc) :: A1, A2 roundoff = 10.*epsilon(roundoff) write(*,*) 'roundoff = ',roundoff !--------------------------------------------------- ! set the arc beg,end geoditic coordinates !--------------------------------------------------- A1%sPnt%lon = 20 ; A1%sPnt%lat = 10 A1%ePnt%lon = 90 ; A1%ePnt%lat = 60 A2%sPnt%lon = 10 ; A2%sPnt%lat = 50 A2%ePnt%lon = 80 ; A2%ePnt%lat = 5 !--------------------------------------------------- ! transfor geoditic coord -> cartesian !--------------------------------------------------- call geod2cart( A1%sPnt%lon, A1%sPnt%lat, & A1%sxyz%x, A1%sxyz%y, A1%sxyz%z ) write(*,*) ' ' write(*,*) 'Start point(lon,lat) arc 1 = ',A1%sPnt%lon,A1%sPnt%lat write(*,*) 'Start point(x,y,z) = ',A1%sxyz%x,A1%sxyz%y, A1%sxyz%z call geod2cart( A1%ePnt%lon, A1%ePnt%lat, & A1%exyz%x, A1%exyz%y, A1%exyz%z ) write(*,*) ' ' write(*,*) 'End point(lon,lat) arc 1 = ',A1%ePnt%lon,A1%ePnt%lat write(*,*) 'End point(x,y,z) = ',A1%exyz%x,A1%exyz%y, A1%exyz%z call geod2cart( A2%sPnt%lon, A2%sPnt%lat, & A2%sxyz%x, A2%sxyz%y, A2%sxyz%z ) write(*,*) ' ' write(*,*) 'Start point(lon,lat) arc 2 = ',A2%sPnt%lon,A2%sPnt%lat write(*,*) 'Start point(x,y,z) = ',A2%sxyz%x,A2%sxyz%y, A2%sxyz%z call geod2cart( A2%ePnt%lon, A2%ePnt%lat, & A2%exyz%x, A2%exyz%y, A2%exyz%z ) write(*,*) ' ' write(*,*) 'End point(lon,lat) arc 2 = ',A2%ePnt%lon,A2%ePnt%lat write(*,*) 'End point(x,y,z) = ',A2%exyz%x,A2%exyz%y, A2%exyz%z call cross_prod( (/A1%sxyz%x,A1%sxyz%y,A1%sxyz%z/), & (/A1%exyz%x,A1%exyz%y,A1%exyz%z/), & N1 ) call cross_prod( (/A2%sxyz%x,A2%sxyz%y,A2%sxyz%z/), & (/A2%exyz%x,A2%exyz%y,A2%exyz%z/), & N2 ) write(*,*) ' ' write(*,*) 'N1 = ',N1(:) write(*,*) 'N2 = ',N2(:) call cross_prod( N1, N2, L ) write(*,*) ' ' write(*,*) 'L = ',L(:) mag = sqrt( sum(L(:)*L(:)) ) I1(:) = L(:)/mag write(*,*) ' ' write(*,*) 'mag = ',mag write(*,*) 'I1 = ',I1(:) radius = sqrt( sum(I1(:)*I1(:)) ) write(*,*) ' ' write(*,*) 'radius = ',radius !--------------------------------------------------- ! Angle test for arc 1, xsect pnt 1 !--------------------------------------------------- write(*,*) ' ' write(*,*) 'Intersect point 1' Spnt(:) = (/A1%sxyz%x,A1%sxyz%y,A1%sxyz%z/) Epnt(:) = (/A1%exyz%x,A1%exyz%y,A1%exyz%z/) Angles(1) = dot_prod_ang( Spnt,I1 ) Angles(2) = dot_prod_ang( Epnt,I1 ) Angles(3) = dot_prod_ang( Spnt,Epnt ) write(*,*) ' ' write(*,*) 'Angles for arc1 = ',Angles(:) tmpsum = Angles(1) + Angles(2) tmp = Angles(3) - tmpsum write(*,*) ' ' write(*,*) 'Check for arc 1 = ',tmpsum,Angles(3), tmp !--------------------------------------------------- ! Angle test for arc 2, xsect pnt 1 !--------------------------------------------------- Spnt(:) = (/A2%sxyz%x,A2%sxyz%y,A2%sxyz%z/) Epnt(:) = (/A2%exyz%x,A2%exyz%y,A2%exyz%z/) Angles(1) = dot_prod_ang( Spnt,I1 ) Angles(2) = dot_prod_ang( Epnt,I1 ) Angles(3) = dot_prod_ang( Spnt,Epnt ) write(*,*) ' ' write(*,*) 'Angles for arc2 = ',Angles(:) tmpsum = Angles(1) + Angles(2) tmp = Angles(3) - tmpsum if( abs(tmp) <= roundoff ) then found = .true. XSECT(:) = I1(:) lambda = asin( I1(3) ) * R2D if( I1(2) /= 0. ) then theta = atan( I1(1)/I1(2) ) * R2D else theta = 90. endif else found = .false. endif write(*,*) ' ' write(*,*) 'Check for arc 1 = ',tmpsum,Angles(3), tmp !--------------------------------------------------- ! Angle test for arc 1, xsect pnt 2 !--------------------------------------------------- if( .not. found ) then I1(:) = -I1(:) write(*,*) ' ' write(*,*) 'Intersect point 2' Spnt(:) = (/A1%sxyz%x,A1%sxyz%y,A1%sxyz%z/) Epnt(:) = (/A1%exyz%x,A1%exyz%y,A1%exyz%z/) Angles(1) = dot_prod_ang( Spnt,I1 ) Angles(2) = dot_prod_ang( Epnt,I1 ) Angles(3) = dot_prod_ang( Spnt,Epnt ) write(*,*) ' ' write(*,*) 'Angles for arc1 = ',Angles(:) !--------------------------------------------------- ! Angle test for arc 2, xsect pnt 2 !--------------------------------------------------- Spnt(:) = (/A2%sxyz%x,A2%sxyz%y,A2%sxyz%z/) Epnt(:) = (/A2%exyz%x,A2%exyz%y,A2%exyz%z/) Angles(1) = dot_prod_ang( Spnt,I1 ) Angles(2) = dot_prod_ang( Epnt,I1 ) Angles(3) = dot_prod_ang( Spnt,Epnt ) write(*,*) ' ' write(*,*) 'Angles for arc2 = ',Angles(:) tmpsum = Angles(1) + Angles(2) tmp = Angles(3) - tmpsum if( abs(tmp) <= roundoff ) then found = .true. XSECT(:) = I1(:) lambda = asin( I1(3) ) * R2D if( I1(2) /= 0. ) then theta = atan( I1(1)/I1(2) ) * R2D else theta = 90. endif endif endif if( found ) then write(*,*) 'Arcs intersect' write(*,*) 'Arcs intersect @ (lon=',theta,',lat=',lambda,')' else write(*,*) 'Arcs do NOT intersect' endif CONTAINS real function dot_prod_ang( V1, V2 ) implicit none real, intent(in) :: V1(3) real, intent(in) :: V2(3) real :: magV1, magV2, angle magV1 = sqrt( sum(V1(:)*V1(:)) ) magV2 = sqrt( sum(V2(:)*V2(:)) ) angle = dot_product( V1,V2 )/(magV1*magV2) dot_prod_ang = acos( angle ) end function dot_prod_ang subroutine cross_prod( A, B, C) implicit none real, intent(in) :: A(3), B(3) real, intent(out) :: C(3) C(1) = A(2)*B(3) - A(3)*B(2) C(2) = A(3)*B(1) - A(1)*B(3) C(3) = A(1)*B(2) - A(2)*B(1) end subroutine subroutine geod2cart( lon, lat, x, y, z ) implicit none real, intent(in) :: lon, lat real, intent(out) :: x, y, z real :: tmp real :: lonr, latr lonr = D2R*lon ; latr = D2R*lat tmp = cos(latr) z = sin(latr) x = tmp*sin(lonr) y = tmp*cos(lonr) end subroutine end program ArcXsect