Re: connect two points by arrow

From: Dave Kennison (kennison@XXXXXX)
Date: Thu Nov 14 2002 - 14:27:07 MST

  • Next message: Laurent Bocahut: "NCAR Graphics on SGI Origin2000"

    The question was as follows:

    > how can I draw an ARROW from point A (alon, alat) to point B (blon,blat)
    > on a map, any suggestion is welcome, thanks.
    >
    > Yanjun

    Appended below is a test program that may do what you want. I have made
    some assumptions that may or may not be correct for your purposes, but I
    think what I have done is the most difficult thing you might want to do:
    it draws the arrow along the shortest great circle route from point A to
    point B, using short segments (the size of the segments being determined
    by a PARAMETER statement in subroutine DRSGCR). The width and length of
    the arrowhead are set in a DATA statement in the main program.

    This was constructed from bits and pieces of code that I had lying about
    and may not be the simplest solution, but I think it does work. Tell me
    if it doesn't do what you want and I will try to modify it accordingly.

    Thanks,

    Dave Kennison

    Program follows:
    ----------------

          PROGRAM TESTIT
    C
    C This program draws an arrow from the point A to the point B on the
    C surface of the globe.
    C
    C NOTE: See the routine DRSGCR for parameters controlling the size of
    C the little segments used to draw the arrow.
    C
    C Define some points A and B to try out.
    C
            DIMENSION ALAT(4),ALON(4),BLAT(4),BLON(4)
    C
            DATA ALAT(1),ALON(1),BLAT(1),BLON(1) / -15., -20., 25., 25. /
            DATA ALAT(2),ALON(2),BLAT(2),BLON(2) / -25., -10., 15., -20. /
            DATA ALAT(3),ALON(3),BLAT(3),BLON(3) / -10., -10., -15., 25. /
            DATA ALAT(4),ALON(4),BLAT(4),BLON(4) / -15., -15., -25., 25. /
    C
    C Define the length and width of the arrowhead.
    C
            DATA ALTH,AWTH / 2. , 1. / ! arrowhead length and width
    C
    C Open GKS.
    C
            CALL OPNGKS
    C
    C Define a color index corresponding to the color red.
    C
            CALL GSCR (1,2,1.,0.,0.)
    C
    C Initialize EZMAP.
    C
            CALL MAPROJ ('OR',0.,0.,0.)
            CALL MAPSET ('AN',40.,40.,40.,40.)
    C
    C Draw a simple map.
    C
            CALL MAPDRW
    C
    C Change the current polyline color to red and double the line width.
    C
            CALL PLOTIT (0.,0.,2)
            CALL GSPLCI (2)
            CALL GSLWSC (2.)
    C
    C Draw the arrows.
    C
            DO 101 I=1,4
              CALL DARROW (ALAT(I),ALON(I),BLAT(I),BLON(I),ALTH,AWTH)
      101 CONTINUE
    C
    C Set the polyline color and line width back to the defaults.
    C
            CALL PLOTIT (0.,0.,2)
            CALL GSPLCI (1)
            CALL GSLWSC (1.)
    C
    C Advance the frame.
    C
            CALL FRAME
    C
    C Close GKS.
    C
            CALL CLSGKS
    C
    C Done.
    C
            STOP
    C
          END

          SUBROUTINE DARROW (ALAT,ALON,BLAT,BLON,ALTH,AWTH)
    C
    C Draw the body of the arrow.
    C
            CALL DRSGCR (ALAT,ALON,BLAT,BLON)
    C
    C Find the latitude and longitude of one end point of the arrowhead
    C and draw it.
    C
            CALL GTLTLN (BLAT,BLON,ALAT,ALON,+AWTH/2.,+ALTH,CLAT,CLON)
            CALL DRSGCR (BLAT,BLON,CLAT,CLON)
    C
    C Find the latitude and longitude of the other end point of the
    C arrowhead and draw it.
    C
            CALL GTLTLN (BLAT,BLON,ALAT,ALON,-AWTH/2.,+ALTH,CLAT,CLON)
            CALL DRSGCR (BLAT,BLON,CLAT,CLON)
    C
    C Done.
    C
            RETURN
    C
          END

          SUBROUTINE GTLTLN (ALAT,ALON,BLAT,BLON,CLAT,CLON,DLAT,DLON)
    C
    C (GTLTLN = GeT LaTitude and LoNgitude)
    C
    C This routine rotates the great circle route between points A and B on
    C the globe so as to put A at the intersection of the prime meridian and
    C the equator and B to the east of A on the equator. It then takes the
    C point at latitude CLAT and longitude CLON in that system and rotates
    C it back to find the point with latitude DLAT and longitude DLON that
    C lies in the same position relative to the original great arc AB as the
    C point with latitude CLAT and longitude CLON does relative to the
    C rotated arc AB.
    C
    C Constants for conversion from degrees to radians and vice-versa.
    C
            DATA DTOR / .017453292519943 /
            DATA RTOD / 57.2957795130823 /
    C
    C Get X, Y, and Z coordinates for point A.
    C
            XCOA=COS(DTOR*ALAT)*COS(DTOR*ALON)
            YCOA=COS(DTOR*ALAT)*SIN(DTOR*ALON)
            ZCOA=SIN(DTOR*ALAT)
    C
    C Get X, Y, and Z coordinates for point B.
    C
            XCOB=COS(DTOR*BLAT)*COS(DTOR*BLON)
            YCOB=COS(DTOR*BLAT)*SIN(DTOR*BLON)
            ZCOB=SIN(DTOR*BLAT)
    C
    C Get X, Y, and Z coordinates for point C.
    C
            XCOC=COS(DTOR*CLAT)*COS(DTOR*CLON)
            YCOC=COS(DTOR*CLAT)*SIN(DTOR*CLON)
            ZCOC=SIN(DTOR*CLAT)
    C
    C Rotate points A and B about the Z axis so as to put point A at
    C longitude zero.
    C
            CALL NGRITD (3,-ALON,XCOA,YCOA,ZCOA)
            CALL NGRITD (3,-ALON,XCOB,YCOB,ZCOB)
    C
    C Rotate points A and B about the Y axis so as to put point A at
    C latitude zero and longitude zero.
    C
            CALL NGRITD (2,+ALAT,XCOA,YCOA,ZCOA)
            CALL NGRITD (2,+ALAT,XCOB,YCOB,ZCOB)
    C
    C Rotate points A and B about the X axis so as to put point B on
    C the equator.
    C
            ANGL=RTOD*ATAN2(ZCOB,YCOB)
    C
            CALL NGRITD (1,-ANGL,XCOA,YCOA,ZCOA)
            CALL NGRITD (1,-ANGL,XCOB,YCOB,ZCOB)
    C
    C Now, inverse rotate the point C.
    C
            CALL NGRITD (1,+ANGL,XCOC,YCOC,ZCOC)
            CALL NGRITD (2,-ALAT,XCOC,YCOC,ZCOC)
            CALL NGRITD (3,+ALON,XCOC,YCOC,ZCOC)
    C
    C Get the latitude and longitude of point D and return them.
    C
            DLAT=RTOD*ASIN(ZCOC)
            DLON=RTOD*ATAN2(YCOC,XCOC)
    C
            RETURN
    C
          END

          SUBROUTINE DRSGCR (ALAT,ALON,BLAT,BLON)
    C
    C (DRSGCR = DRaw Shortest Great Circle Route)
    C
    C This routine draws the shortest great circle route joining two points
    C on the globe.
    C
            PARAMETER (MPTS=360,SIZE=.5) ! MPTS = INT(180./SIZE) + 1
    C
            DIMENSION QLAT(MPTS),QLON(MPTS)
    C
            NPTS=MAX(1,MIN(MPTS,INT(ADSGCR(ALAT,ALON,BLAT,BLON)/SIZE)))
    C
            CALL MAPGCI (ALAT,ALON,BLAT,BLON,NPTS,QLAT,QLON)
    C
            CALL MAPIT (ALAT,ALON,0,IDIS)
    C
            DO 101 I=1,NPTS
              CALL MAPIT (QLAT(I),QLON(I),1)
      101 CONTINUE
    C
            CALL MAPIT (BLAT,BLON,1)
    C
            CALL MAPIQ
    C
            RETURN
    C
          END

          FUNCTION ADSGCR (ALAT,ALON,BLAT,BLON)
    C
    C (ADSGCR = Angle in Degrees along Shortest Great Circle Route)
    C
    C This function returns the shortest great circle distance, in degrees,
    C between two points, A and B, on the surface of the globe.
    C
            DATA DTOR / .017453292519943 /
            DATA RTOD / 57.2957795130823 /
    C
            XCOA=COS(DTOR*ALAT)*COS(DTOR*ALON)
            YCOA=COS(DTOR*ALAT)*SIN(DTOR*ALON)
            ZCOA=SIN(DTOR*ALAT)
    C
            XCOB=COS(DTOR*BLAT)*COS(DTOR*BLON)
            YCOB=COS(DTOR*BLAT)*SIN(DTOR*BLON)
            ZCOB=SIN(DTOR*BLAT)
    C
            ADSGCR=2.*RTOD*ASIN(SQRT((XCOA-XCOB)**2+
         + (YCOA-YCOB)**2+
         + (ZCOA-ZCOB)**2)/2.)
    C
            RETURN
    C
          END

          SUBROUTINE NGRITD (IAXS,ANGL,UCRD,VCRD,WCRD)
    C
    C Define a multiplicative constant to convert from degrees to radians.
    C
            DATA DTOR / .017453292519943 /
    C
    C This routine rotates the point with coordinates (UCRD,VCRD,WCRD) by
    C the angle ANGL about the axis specified by IAXS (1 for the U axis,
    C 2 for the V axis, 3 for the W axis). A right-handed coordinate
    C system is assumed.
    C
            SINA=SIN(DTOR*ANGL)
            COSA=COS(DTOR*ANGL)
    C
            UTMP=UCRD
            VTMP=VCRD
            WTMP=WCRD
    C
            IF (IAXS.EQ.1) THEN
              VCRD=VTMP*COSA-WTMP*SINA
              WCRD=WTMP*COSA+VTMP*SINA
            ELSE IF (IAXS.EQ.2) THEN
              UCRD=UTMP*COSA+WTMP*SINA
              WCRD=WTMP*COSA-UTMP*SINA
            ELSE
              UCRD=UTMP*COSA-VTMP*SINA
              VCRD=VTMP*COSA+UTMP*SINA
            END IF
    C
            RETURN
    C
          END
    _______________________________________________
    ncarg-talk mailing list
    ncarg-talk@ucar.edu
    http://mailman.ucar.edu/mailman/listinfo/ncarg-talk



    This archive was generated by hypermail 2b29 : Thu Nov 14 2002 - 14:37:29 MST