From: [email-masked] (Gonzalo Hernandez)
Date: Mon, 20 Nov 89 22:15:13 PST

C Ben:
C  This program uses the data from rayplot and sets it up for a
C  hairy-plot.  You should comment out all the calls for the plotting,
C  ( i.e., call gpaper , call plabel, call line, etc. ).  Just make sure
C  that the IA parameter is properly set as 0.  The output should look
C  like the test file also sent for your testing your troff.
C
      PROGRAM RAYGRAPS
C *** Usage: link data to be plotted as fort2
C *** Grap output in fort9
C *** Plots wind data as vectors as a function of time on a sun dial.
C *** Southern Hemisphere version....
C *** This program needs the evenly spaced spline data from Ray's 
C *** model in MLT and 25 points.
C *** Statements 1->9; 12; 28; 51. Formats 100->104; 107->109,
C *** 113; 200->201; 203; 210->211; 300 -> 318
      DIMENSION XA(2),YA(2),TIMEUT(12,30),VN(12,30),UN(12,30),
     *UMAG(12,30),GEOMLAT(12),IXP(6)
      COMMON/HAIR/PX(600),Y(600)
      COMMON /GSIZE/ XP1,XP2,YP1,YP2,XNEW,YNEW,SIZELABX,SIZELABY,
     *XLABDISP,YLABDISP,TICKHEIGHT,BIGLABSIZE,BIGLABDISP,MAXWIDTHY,
     *MAXHEIGHTX,IAX,IAXX,IAY,IAYY,ISTATION,IFRA
      COMMON /GTITLE/ ITOP,IX,IY,IXF,IYF
      COMMON /GBOUND/ GX(3),GY(3)
      COMMON /GLABEL/ XLAB(50),YLAB(50),NUMX,NUMY
      CHARACTER*2 ID
      CHARACTER*10 IXF,IYF,ILAB
      CHARACTER*65 ITOP, IX, IY, LLABEL
      CHARACTER*80 LINE(17)
      CHARACTER*3 ISTAT, ISTA, LAB
      CHARACTER*4 XLIN,YLIN,XLOG,YLOG,XGRD,YGRD,XTIK,YTIK
      DATA XLIN /4HXLIN/,XLOG /4HXLOG/,XGRD /4HXGRD/,XTIK /
     & 4HXTIK/, YLIN /4HYLIN/, YLOG /4HYLOG/,YGRD /4HYGRD/,YTIK/4HYTIK/
      DATA PI / 3.141592654 /
      CALL GPSET
      BIGLABSIZE = 3.
      WRITE(*,102)
  102 FORMAT(1X,' Hairy wind line plotting program ',/ )
      WRITE(*,200)
  200 FORMAT(1X,' Type the lowest |latitude| desired (degrees): ')
      READ(*,*) RADIUS
      GMLAT = RADIUS
      RADIUS = 90. - RADIUS
      WRITE(*,203)
  203 FORMAT(1X,' Choose plotting station ( 0, 1, 3): ')
      READ(*,*) ISTATION
      BSIZE = 4.0
      WRITE(*,300)
  300 FORMAT(1X,'Please type the box side dimension in inches: ')
      READ(*,*)BSIZE
    1 J = 0
    6 READ(2,100,END=7)IADATE,IDAY,IYR,ISTAT,GEOMLAT(J+1),
     *(TIMEUT(J+1,I),UMAG(J+1,I),VN(J+1,I),UN(J+1,I), I = 1,25 )
  100 FORMAT(1X,3I10,A3,F10.3,/,(2F7.2,2F10.2))
      IF(GEOMLAT(J+1) .LT. GMLAT) GO TO 6
      J = J + 1
      IF(J .LT. 12) GO TO 6
    7 CONTINUE
      NPJ = J
      IF( NPJ .LT. 1 ) GO TO 5
      NPI = 24
      WRITE(ISTA,201)ISTAT
  201 FORMAT(A3)
  101 FORMAT(1X,F10.2,5F10.2,F10.3)
C *** Find the longest vector
      TOPM = -1.0E30
      TOPZ = TOPM
      TOPVT = TOPM
      DO 4 J = 1, NPJ
      DO 4 I = 1, NPI
      TOPM = AMAX1(TOPM,ABS(VN(J,I)))
    4 CONTINUE
      SCALEW = 100.
      TOPT = AINT(TOPM/SCALEW)*SCALEW + SCALEW
      WRITE(*,104) TOPT
  104 FORMAT(1X,'Largest wind is ',F5.0,' Type desired scale ',
     *'(-9 to ignore ): ')
      READ(*,*)WVMAX
      IF(WVMAX .LT. -1.0) GO TO 9
      TOPT = WVMAX
    9 CONTINUE
      TOPVT = TOPT
      ITOPT = INT(TOPT)
      WSCALE = 3.5/TOPT
      RADIUS = RADIUS + 3.5
   12 WRITE(IXF,210)
  210 FORMAT(10H(F3.0,1H#))
      WRITE(IYF,211)
  211 FORMAT(10H(F5.0,1H#))
      YLAB( 1 ) = -9999.
      XLAB( 1 ) = -9999.
      GY(2)= RADIUS
      GY(1)= -GY(2)
      GY(3)= GY(2)
      GX(2) = GY(2)
      GX(1) = -GY(2)
      GX(3) = GY(2)
C*****grap data
      WRITE(LINE(1),301)
  301 FORMAT(3H.G1,77X)
C *** Box size 4-liner
      BRAD = BSIZE/(2.0*RADIUS)
      WRITE(LINE(2),302)BSIZE,BSIZE
  302 FORMAT(9Hframe ht ,F6.2,5H wid ,F6.2, 54X)
      WRITE(LINE(3),303) GX(1),GX(2),GY(1),GY(2)
  303 FORMAT(8Hcoord x ,F7.3,2H ,,F7.3,3H y ,F7.3,2H ,,F7.3,35X)
C     WRITE(LINE(4),304)BRADIUS
  304 FORMAT(22Hcircle at 0, 0 radius ,F7.3,51X)
C     WRITE(LINE(5),305)BRADIUS
  305 FORMAT(22Hcircle at 0, 0 radius ,F7.3,51X)
      WRITE(LINE(6),306)ISTA,IADATE,IDAY
  306 FORMAT(10Hlabel top ,4H"\fB,A3,8H 15867K ,5X,I6,1X,
     *4H(Day,I4,5H)\fR" ,8H size +2 , 22X)
      XA(1) = 0.0
      YA(1) = -RADIUS*0.98
      WRITE(LINE(7),307)XA(1),YA(1)
  307 FORMAT(19H"\fB0\fR" above at , F7.2,1H,,F7.2, 46X)
      XA(1) = -RADIUS*0.9245
      YA(1) =  0.0
      WRITE(LINE(8),308)XA(1),YA(1)
  308 FORMAT(19H"\fB6\fR" rjust at , F7.2,1H,,F7.2, 46X)
      XA(1) = 0.0
      YA(1) = RADIUS*0.9608
      WRITE(LINE(9),309)XA(1),YA(1)
  309 FORMAT(20H"\fB12\fR" below at , F7.2,1H,,F7.2, 45X)
      XA(1) =  RADIUS*0.891
      YA(1) =  0.0
      WRITE(LINE(10),310)XA(1),YA(1)
  310 FORMAT(20H"\fB18\fR" ljust at , F7.2,1H,,F7.2, 45X)
C     WRITE(LINE(11),311)KLABEL,ZERO,RSTARTN
  311 FORMAT(4H"\fB,I2,15H S\fR" above at ,F7.2,1H,,F7.2,43X)
C     WRITE(LINE(12),312)XA(1),YA(1),XA(2),YA(2)
  312 FORMAT(11Harrow from ,F7.2,2H, ,F7.2,4H to ,F7.2,2H, ,F7.2,35X)
      WRITE(LINE(13),313)
  313 FORMAT(3H.G2 ,77X)
      WRITE(LINE(14),314)
  314 FORMAT(3H.bp,77X)
      WRITE(9,315)(LINE(I),I=1,3)
      WRITE(LINE(16),316)
  316 FORMAT(14Hticks left off, 66X)
      WRITE(LINE(17),317)
  317 FORMAT(16Hticks bottom off, 64X)
      WRITE(9,315) LINE(16), LINE(17)
      WRITE(9,315)(LINE(I),I=6,10)
  315 FORMAT(A80)
      XA(1) = RADIUS - 7.0
      XA(2) = XA(1) + TOPT*WSCALE
      YA(1) = -RADIUS + 3.5
      YA(2) = -RADIUS + 3.5
      WRITE(LINE(12),312)XA(1),YA(1),XA(2),YA(2)
      XA(1) = XA(1) + TOPT*WSCALE/2.0
      WRITE(LINE(15),318) ITOPT,XA(1),YA(2)
  318 FORMAT(1H",I4,23H m/s" size -3 above at ,F8.2,1H,,F8.2, 36X )
      WRITE(9,315) LINE(12)
      WRITE(9,315) LINE(15)
      WRITE(ITOP,113)ISTA,IADATE,IDAY
  113 FORMAT(17X,A3,1X, 7H15867K ,4X,I6,1X,4H(Day,I4,1H),16X)
      CALL GPAPER(XLIN,YLIN,XTIK,YTIK,IA)
      IF(IA .LT. 0 ) GO TO 28
      PCY = 50. - 0.5
      PCX = 1.75/RADIUS*50. -0.5
      KLABEL = 6
      WRITE(LLABEL,107)KLABEL
  107 FORMAT(31X,I2,31X)
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      PCX = (2.0*RADIUS -1.75)/RADIUS*50. - 2.0
      KLABEL = 18
      WRITE(LLABEL,107)KLABEL
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      PCX = 50. - 0.9
      PCY = 1.75/RADIUS*50. - 0.5
      KLABEL = 0
      WRITE(LLABEL,107)KLABEL
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      PCX = 50. - 1.25
      PCY = (2.0*RADIUS -1.75)/RADIUS*50. - 1.0
      KLABEL = 12
      WRITE(LLABEL,107)KLABEL
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      RSTART = 11.5
      DX = RSTART/50.
      CALL HAIRCIR (RSTART, DX, NPT )
      CALL LINES(PX,Y,NPT,3,IA)
      IF(IA .LT. 0 ) GO TO 28
      BRADIUS = BRAD*RSTART
      WRITE(LINE(5),305)BRADIUS
C     WRITE(9,315) LINE(5)
      RSTART = 15.0
      DX = RSTART/50.
      CALL HAIRCIR (RSTART, DX, NPT )
      CALL LINES(PX,Y,NPT,3,IA)
      IF(IA .LT. 0 ) GO TO 28
      BRADIUS = BRAD*RSTART
      WRITE(LINE(5),305)BRADIUS
C     WRITE(9,315) LINE(5)
      RSTART = 18.5
      DX = RSTART/50.
      CALL HAIRCIR (RSTART, DX, NPT )
      CALL LINES(PX,Y,NPT,3,IA)
      IF(IA .LT. 0 ) GO TO 28
      BRADIUS = BRAD*RSTART
      WRITE(LINE(5),305)BRADIUS
C     WRITE(9,315) LINE(5)
      RSTART = 5.0
    3 DX = RSTART/50.
      CALL HAIRCIR (RSTART, DX, NPT )
      CALL LINES(PX,Y,NPT,3,IA)
      IF(IA .LT. 0 ) GO TO 28
      BRADIUS = BRAD*RSTART
      WRITE(LINE(4),304)BRADIUS
      WRITE(9,315) LINE(4)
      PSIZE = 1.
      PCX = 50. - 1.25
      PCY = (RADIUS -RSTART)/RADIUS*50. - 0.5
      KLABEL = INT(90. - RSTART)
      WRITE(LLABEL,107)KLABEL
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      RSTARTN = - RSTART
      ZERO = 0.0
      WRITE(LINE(11),311)KLABEL,ZERO,RSTARTN
      WRITE(9,315) LINE(11)
      RSTART = RSTART + 5.0
      IF( RSTART .LT. RADIUS ) GO TO 3
      ITOPVT = INT(TOPVT)
      WRITE(LLABEL,109)ITOPVT
  109 FORMAT(28X,I3,4H m/s,29X)
      PSIZE = 1.
      PCX = 79.
      PCY = 92.
      CALL PLABEL(PCX,PCY,PSIZE,LLABEL,IH,IV,IA)
      IF(IA .LT. 0 ) GO TO 28
      YA(1) = 0.80*RADIUS
      YA(2) = YA(1)
      VT = TOPVT*WSCALE/2.0
      XA(1) = 0.68*RADIUS - VT
      XA(2) = 0.68*RADIUS + VT
      CALL LINES(XA,YA,2,1,IA)
      IF(IA .LT. 0 ) GO TO 28
      DO 8 J = 1, NPJ
        DO 2 I = 1, NPI
        VN(J,I) = VN(J,I)*WSCALE
        UN(J,I) = UN(J,I)*WSCALE
    2   CONTINUE
      RSTART = 90. - GEOMLAT(J)
       CALL HAIRLINE(VN, UN, UMAG, RSTART, NPI, J )
    8 CONTINUE
      WRITE(9,315) LINE(13)
      WRITE(9,315) LINE(14)
C     GO TO 1
      GO TO 5
   28 WRITE(*,103) IA
  103 FORMAT(1X, 'IA .LT. 0, i.e, ',I4,' PROGRAM END ')
      GO TO 5
   51 WRITE(*,108)
  108 FORMAT(1X,'READING OUT-OF-SEQUENCE DATA')
    5 CONTINUE
      CALL GPEND 
      END
      SUBROUTINE HAIRCIR (R, DX, NPT )
      COMMON/HAIR/PX(600),Y(600)
      DATA PI / 3.141592654 /
      NPI = INT(R/DX)
      NPC = NPI + 1
      NP2 = 2*NPI + 2
      NP3 = 2*NPI
      NP4 = 4*NPI + 2
      ANGLE = PI/(2.0*NPI)
      DO 2 I = 1, NPC
      ANGLES = ANGLE*FLOAT(I - 1)
      X  = R * SIN(ANGLES)
      PX(I) = X
      Y(I)  = R * COS(ANGLES)
      PX(NP2 - I) = X
      Y( NP2 - I) = -Y(I)
      PX(NP3 + I) = -X
      Y( NP3 + I) = -Y(I)
      PX(NP4 - I) = -X
      Y( NP4 - I) = Y(I)
    2 XX = XX + DXX
      NPT = 4*NPI + 1
      RETURN
      END
      SUBROUTINE HAIRLINE(VM, VZ, TIME, R, NP, J )
      DIMENSION VM(12,1), VZ(12,1), TIME(12,1), X(2), Y(2)
      CHARACTER*80 LINE(17)
      DATA PI180 / 0.017453293 /
      DO 1 I = 1, NP
      ANGLE = ( 270. - (TIME(J,I)*15) )*PI180
      X(1) = COS(ANGLE)
      Y(1) = SIN(ANGLE)
      X(2) = (R + VM(J,I))*X(1) + VZ(J,I)*Y(1)
      Y(2) = (R + VM(J,I))*Y(1) - VZ(J,I)*X(1)
      X(1) = X(1)*R
      Y(1) = Y(1)*R
      LINETYPE = 1
      CALL LINES( X,Y,2,LINETYPE,IA)
      IF(IA .LT. 0 ) GO TO 1
      WRITE(LINE(12),312)X(1),Y(1),X(2),Y(2)
  312 FORMAT(18Harrow dashed from ,F7.2,2H, ,F7.2,4H to ,F7.2,2H, ,F7.2,
     *27X)
      WRITE(9,315) LINE(12)
  315 FORMAT(A80)
    1 CONTINUE
      RETURN
      END