;+
; Project     : Course Homework Program.
;
; Name        : PROJ_TUTORIAL
;
; Purpose     : Procedure to calculate projectiles by a variety of methods.
;
; Explanation : This procedure is used as a tutorial to calculate 
;               projectiles analytically by algebraic equations and by the
;               solution of ODEs using both the Euler Method and RK4.  It 
;               asks the user to enter an initial velocity and projection 
;               angle of a projectile.  This code also has the option of 
;               inputting the initial time, position, and height of launch.
;               Results are compared with each method of solution.  Plots  
;               are made of each solution.
;
; Use         : PROJ_TUTORIAL (no passed parameters).
;
; Inputs      : None.
;
; Opt. Inputs : None.
;
; Outputs     : None.
;
; Opt. Outputs: None.
;
; Keywords    : VINIT:  The inital velocity in meters/sec.
;
;               ANGINIT:  The projection angle of launch in degrees.
;
;               TINIT:  The time of launch.  [DEFAULT: 0. sec]
;
;               XINIT:  The initial horizontal position.  [DEFAULT:  0. m]
;
;               YINIT:  The initial vertical position.  [DEFAULT:  0. m]
;
;               TSTEP:  The time step used in the solution to the ODEs.
;                       [DEFAULT:  0.1 s]
;
;               GRAV:   The surface gravity.  [DEFAULT:  9.8038 m/s^2]
;
;               ERRMSG: If defined and passed, then any error messages will
;                       be returned to the user in this parameter rather than
;                       being printed internally by this procedure.  If no
;                       errors are encountered, then a null string is
;                       returned.  In order to use this feature, the string
;                       ERRMSG must be defined first, e.g.,
;
;                               ERRMSG = ''
;                               PROJ_TUTORIAL, ERRMSG=ERRMSG
;                               IF ERRMSG[0] NE '' THEN ...
;
; Calls       : DIALOG_PROJ_INPUT, PROJ_RK4_EQN.
;
; Common      : GRAVSURF.
;
; Restrictions: None.
;
; Side effects: None.
;
; Category    : Class tutorial.
;
; Prev. Hist. : None.
;
; Written     : Donald G. Luttermoser, ETSU/Physics, 3 November 2015.
;
; Modified    : Version 1, Donald G. Luttermoser, ETSU/Physics, 3 Nov 2015
;                       Initial program.
;
; Version     : Version 1,  3 November 2015.
;
;-
PRO PROJ_TUTORIAL, ERRMSG=ERRMSG
;
COMMON GRAVSURF, GRAV
;
ON_ERROR, 2     ; Return to caller if error is encountered.
EMESS = ''      ; Set to non-null string if error is encountered.
;
; Define generic internal parameters.
;
PASS = 0      ; PASS = 0: Plot to terminal;  PASS > 0: to hardcopy.
AFONT = [-1, 0]
TXTSIZE = [1.2, 1.5]
ASK = 'Y'
STHETA = ['!7h', '!7q']
TERM = STRUPCASE(!D.NAME)
PLNAMES = ['projectile_pos.ps', 'projectile_yt.ps', 'projectile_vyt.ps'] 
;
; Check for keywords, force the user to use double precision to cut 
; down on round-off error.
;
;VINIT, ANGINIT, TINIT, XINIT, YINIT, TSTEP, GRAV
IF KEYWORD_SET(VINIT) THEN V0 = DOUBLE(VINIT) ELSE V0 = 0.D0
IF KEYWORD_SET(ANGINIT) THEN PANGLE = DOUBLE(ANGINIT) ELSE PANGLE = 0.D0
IF KEYWORD_SET(TINIT) THEN T0 = DOUBLE(TINIT) ELSE T0 = 0.D0
IF KEYWORD_SET(XINIT) THEN X0 = DOUBLE(XINIT) ELSE X0 = 0.D0
IF KEYWORD_SET(YINIT) THEN Y0 = DOUBLE(YINIT) ELSE Y0 = 0.D0
IF NOT KEYWORD_SET(TSTEP) THEN TSTEP = 0.1D0
IF NOT KEYWORD_SET(GRAV)  THEN GRAV = 9.8038D0
;
; User input of starting values.
;
QINPUT = {INPUT, V0: V0, PANGLE: PANGLE, T0: T0, X0: X0, Y0: Y0, $
    TSTEP: TSTEP, GRAV: GRAV}
RESULT = DIALOG_PROJ_INPUT(QINPUT)
IF RESULT EQ 1 THEN BEGIN
    EMESS = 'User stopped program on data input.'
    GOTO, HANDLE_ERROR
ENDIF
V0 = QINPUT.V0  &  PANGLE = QINPUT.PANGLE
T0 = QINPUT.T0  &  TSTEP = QINPUT.TSTEP
X0 = QINPUT.X0  &  Y0 = QINPUT.Y0
GRAV = QINPUT.GRAV
;
; Algebraic solutions:
;
; Calculate the initial vector components (note that !DTOR) converts
; degrees to radians):
;
V0X = V0 * COS(PANGLE * !DTOR)
V0Y = V0 * SIN(PANGLE * !DTOR)
;
; Calculate the maximum height using vy^2 = v0y^2 - 2 g (y - y0) and
; noting that vy = 0 at the top of the trajectory:
;
YMAX = Y0 + (V0Y^2 / (2 * GRAV))
;
; Calculate time to get to YMAX, the total time is then just double
; that time since the problem is symmetrical about that point since
; we are ignoring air friction.  To get time to YMAX, use
; vy = v0y - g * (t - t0), remembering that vy = 0 at the top.
;
T2YMAX = T0 + (V0Y / GRAV)
T_TOTAL = 2 * T2YMAX
;
; XMAX now just results from x = x0 + v0x * (t - t0):
;
XMAX = X0 + V0X * (T_TOTAL - T0)
;
; Now print the required output to the terminal.
;
V0STR = STRTRIM(STRING(V0, FORMAT='(F12.4)'), 2)
A0STR = STRTRIM(STRING(PANGLE, FORMAT='(F8.4)'), 2)
YMAXSTR = STRTRIM(STRING(YMAX, FORMAT='(F12.4)'), 2)
XMAXSTR = STRTRIM(STRING(XMAX, FORMAT='(F12.4)'), 2)
TTOTSTR = STRTRIM(STRING(T_TOTAL, FORMAT='(F12.4)'), 2)
PRINT, 'ALGEBRAIC SOLUTION:'
PRINT, 'For an initial velocity of ', V0STR, ' m/s'
PRINT, ' and an initial projection angle of ', A0STR, ' degrees,'
PRINT, ' '
PRINT, ' this projectile will reach a height of ', YMAXSTR, ' meters,'
PRINT, ' and travel a distance downrange of ', XMAXSTR, ' meters'
PRINT, ' in a time of ', TTOTSTR, ' seconds.'
;
XPATH = FINDGEN(1001) * XMAX / 1000.
TPATH = FINDGEN(1001) * T_TOTAL / 1000.
YPATH = Y0 + V0Y * (TPATH - T0) - 0.5 * GRAV * (TPATH - T0)^2
VYPATH = V0Y - GRAV * (TPATH - T0)
;
; Calculate the trajectory from the ODEs using the Euler method.
;
; Note that our 4 differential equations (DYDT) are
;     DYDT(0) = dx/dt  = vx
;     DYDT(1) = dvx/dt = 0
;     DYDT(2) = dy/dt  = vy
;     DYDT(3) = dvy/dt = -g
;
; Now integrate using the Euler method:
;
TCURR = T0  &  YCURR = Y0  &  VXCURR = V0X  &  VYCURR = V0Y
TEULER = T0  &  XEULER = X0  &  YEULER = Y0
VXEULER = V0X  &  VYEULER = V0Y
XLAST = X0  &  YLAST = Y0  &  TLAST = T0
WHILE TCURR LE T_TOTAL DO BEGIN
    XCURR = XLAST + VXCURR*TSTEP
    YCURR = YLAST + VYCURR*TSTEP
    VYCURR = VYCURR - GRAV*TSTEP
    TCURR = TLAST + TSTEP
    TEULER = [TEULER, TCURR]
    XEULER = [XEULER, XCURR]  &  YEULER = [YEULER, YCURR]
    VXEULER = [VXEULER, VXCURR]  &  VYEULER = [VYEULER, VYCURR]
    TLAST = TCURR  &  XLAST = XCURR  &  YLAST = YCURR
ENDWHILE
;
; Note that our 4 corresponding dependent variables (Y) are:
;     Y(1) = x(i),   Y(2) = vx(i)
;     Y(3) = y(i),   Y(4) = vy(i)
;
; Now integrate using the 4th-order Runge-Kutta method.
; 
POSRK = DBLARR(2)  &  VELRK = POSRK ; Position & velocity in x & y coordinates for RK4.
POSRK[0] = X0  &  POSRK[1] = Y0     ; x & y positions have units of meters.
VELRK[0] = V0X     &  VELRK[1] = V0Y  ; Velocities have units of m/s.
TCURR = T0  &  NTOT = FIX((T_TOTAL - T0) / TSTEP) + 1
; Create arrays for RK4 calculations.
TRK = DBLARR(NTOT)  &  XRK = TRK  &  YRK = TRK  &  VXRK = TRK  &  VYRK = TRK
; Set counter and initial array values.
ISTEP = 1  &  TRK[0] = T0  &  XRK[0] = X0  &  YRK[0] = Y0
VXRK[0] = V0X  &  VYRK[0] = V0Y
WHILE ISTEP LT NTOT DO BEGIN
    STATE = [POSRK, VELRK]
    TCURR = TCURR + TSTEP
    DYDT = PROJ_RK4_EQN(TCURR, STATE, ERRMSG=EMESS)
    IF EMESS NE '' THEN GOTO, HANDLE_ERROR
    RESULT = RK4(STATE, DYDT, TRK(ISTEP), TSTEP, 'proj_rk4_eqn', /DOUBLE)
    POSRK[0] = RESULT[0]  &  POSRK[1] = RESULT[1]
    VELRK[0] = RESULT[2]  &  VELRK[1] = RESULT[3]
;
; Store data in output/plotting arrays.
;
    XRK[ISTEP] = POSRK[0]   &  YRK[ISTEP] = POSRK[1]  &  TRK[ISTEP] = TCURR
    VXRK[ISTEP] = VELRK[0]  &  VYRK[ISTEP] = VELRK[1]
    ISTEP = ISTEP + 1
ENDWHILE
;
; Select the type of plot you want.
;
PLTYPE = ['y(x) plot', 'y(t) plot', 'vy(t) plot']
QPLOT = DIALOG_SELECT(PLTYPE, TITLE='Plot Type', /EXCLUSIVE)
CASE QPLOT OF
    PLTYPE[1]:  NTPL = 1
    PLTYPE[2]:  NTPL = 2
    ELSE:       NTPL = 0
ENDCASE
;
; Plot variables:
;
CASE NTPL OF
    1:  BEGIN  ; y(t)
        XPLOT = TPATH  &  YPLOT = YPATH
        XEPLOT = TEULER  &  YEPLOT = YEULER
        XRKPLOT = TRK  &  YRKPLOT = YRK
        XTIT = '!6Time (s)'  &  YTIT = '!6Height (m)'
        END
    2:  BEGIN  ; vy(t)
        XPLOT = TPATH  &  YPLOT = VYPATH
        XEPLOT = TEULER  &  YEPLOT = VYEULER
        XRKPLOT = TRK  &  YRKPLOT = VYRK
        XTIT = '!6Time (s)'  &  YTIT = '!6Vertical Velocity (m/s)'
        END
    ELSE:  BEGIN  ; y(x)
        XPLOT = XPATH  &  YPLOT = YPATH
        XEPLOT = XEULER  &  YEPLOT = YEULER
        XRKPLOT = XRK  &  YRKPLOT = YRK
        XTIT = '!6Horizontal Distance (m)'  &  YTIT = '!6Height (m)'
        END
ENDCASE
;
; Y-axis range.
;
YMIN = MIN([YPLOT, YEPLOT, YRKPLOT])  &  YMAX = MAX([YPLOT, YEPLOT, YRKPLOT])
YDIFF = YMAX - YMIN
;
; Set box characteristics
;
!X.THICK = 2.4  &  !Y.THICK = 2.4  &  !P.THICK = 2.4
!P.CHARSIZE = 1.5  &  !X.CHARSIZE = 1.0  &  !Y.CHARSIZE = 1.0
;
PLOTIT:
;
PLOT, XPLOT, YPLOT, XTITLE=XTIT, YTITLE=YTIT, FONT=AFONT[PASS], $
    CHARSIZE=TXTSIZE[PASS], THICK=2*(2*PASS+1), YRANGE=[YMIN, YMAX], $
    TITLE='!17Projectile Motion', XCHARSIZE=1.0, YCHARSIZE=1.0
OPLOT, XEPLOT, YEPLOT, THICK=2*(2*PASS+1), LINE=2
OPLOT, XRKPLOT, YRKPLOT, LINE=3, PSYM=2
XYOUTS, 0.5*MAX(XPLOT), YMIN+0.3*YDIFF, '!6v!Do!N = '+V0STR+' m/s', $
    ALIGN=0.5, SIZE=1.4, FONT=AFONT[PASS]
XYOUTS, 0.5*MAX(XPLOT), YMIN+0.2*YDIFF, STHETA[PASS]+'!6!Do!N = '+ $
    A0STR+'!Uo!N', ALIGN=0.5, SIZE=1.4, FONT=AFONT[PASS]
;
XYOUTS, 0.5*MAX(XPLOT), YMIN+0.7*YDIFF, '!6Solid: Analytic', ALIGN=0.5, $
    SIZE=1.4, FONT=AFONT[PASS]
XYOUTS, 0.5*MAX(XPLOT), YMIN+0.6*YDIFF, '!6Dash: Euler Method', ALIGN=0.5, $
    SIZE=1.4, FONT=AFONT[PASS]
XYOUTS, 0.5*MAX(XPLOT), YMIN+0.5*YDIFF, '!6Star: RK4 Method', ALIGN=0.5, $
    SIZE=1.4, FONT=AFONT[PASS]
;
IF PASS EQ 0 THEN BEGIN
    ASK = DIALOG_MESSAGE('Save plot as hardcopy?', /QUESTION)
    IF STRUPCASE(ASK) EQ 'YES' THEN ASK = 'Y' ELSE ASK = 'N'  
ENDIF ELSE BEGIN
;
; Close postscript file.
;
    PRINT, ' '
    PRINT, 'Your plot is now stored in file: ' + PLNAMES[NTPL]
    PRINT, ' '
    DEVICE, /CLOSE
    ASK = 'N'
ENDELSE
;
IF STRUPCASE(ASK) NE 'N' THEN BEGIN
;
; Set the postscript fonts.
;
    PASS = PASS + 1
    SET_PLOT, 'PS'
    DEVICE, /LANDSCAPE, FILENAME=PLNAMES[NTPL]
    DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
    DEVICE, /TIMES, /BOLD, FONT_INDEX=6
    DEVICE, /COURIER, /BOLD, FONT_INDEX=5
    DEVICE, /SYMBOL, FONT_INDEX=7
;
	IF PASS EQ 1 THEN GOTO, PLOTIT
ENDIF
;
SET_PLOT, TERM
;
RETURN
;
; Error handling portion of the procedure.
;
HANDLE_ERROR:
;
IF N_ELEMENTS(ERRMSG) EQ 0 THEN PRINT, EMESS
ERRMSG = EMESS
RETURN
;
END