c
c=====================================================================
c
	program choreograph
c
c=====================================================================
c
	implicit none
c
c----- Initial camera and renderer information (at start of shot)
c
	real xfromi, yfromi, zfromi
	real xati, yati, zati
	real xfi, yfi, zfi
	real fnorm
	real xui, yui, zui
c
c----- Current camera and renderer information (at start of frame)
c
	real xfromc, yfromc, zfromc
	real xatc, yatc, zatc
	real xfc, yfc, zfc
	real xuc, yuc, zuc
c
c----- Frame camera and renderer information (as modified by calculation)
c
	real xfromf, yfromf, zfromf
	real xatf, yatf, zatf
	real xff, yff, zff
	real xuf, yuf, zuf
c	
c----- Translation of from point
c
	logical trans_from, trans_from_fc, trans_from_pc
	real trans_from_vec_x, trans_from_vec_y, trans_from_vec_z
	real trans_from_pt_x, trans_from_pt_y, trans_from_pt_z
	real trans_from_dist_per_frame, trans_from_dist
	real trans_from_dx, trans_from_dy, trans_from_dz
c	
c----- Translation of at point
c
	logical trans_at, trans_at_fc
	real trans_at_vec_x, trans_at_vec_y, trans_at_vec_z
	real trans_at_dist_per_frame, trans_at_dist
	real trans_at_dx, trans_at_dy, trans_at_dz
c
c----- Circle maneuver
c
	logical circle, circle_atc, circle_uc
	real circle_pt_x, circle_pt_y, circle_pt_z
	real circle_vec_x, circle_vec_y, circle_vec_z
	real circle_degrees, circle_degrees_per_frame
c
c----- Pan maneuver
c
	logical pan, pan_uc
	real pan_degrees, pan_degrees_per_frame
	real pan_vec_x, pan_vec_y, pan_vec_z
c
c----- Tilt maneuver
c
	logical tilt, tilt_uc, tilt_sc
	real tilt_degrees, tilt_degrees_per_frame
	real tilt_vec_x, tilt_vec_y, tilt_vec_z
c
c----- Twist maneuver
c
	logical twist, twist_uc, twist_fc
	real twist_degrees, twist_degrees_per_frame
	real twist_vec_x, twist_vec_y, twist_vec_z
c
c----- Other
c
	integer ichoice, num, nframes
	real xtemp, ytemp, ztemp
	real x1, y1, z1
	real x2, y2, z2
	real zmult
	character ans*1
c
	integer iframe, jf
	character outfname*40
c
c=====================================================================
c
c----- GET INITIAL CAMERA AND/OR RENDERER INFORMATION
c
	print*,' Enter initial camera and render information'
c
c----- Get initial frame number
c
	print*,' -- Enter initial frame number (normally = 0)'
	read*, iframe
c
c----- Get from point (camera position)
c
	print*,' -- Enter from point (camera position)'
	read*, xfromi, yfromi, zfromi
c
c----- Get at point (camera look-at point)
c
	print*,' -- Enter at point (camera look-at point)'
	read*, xati, yati, zati
c
c----- Calculate facing vector (vector from to at)
c
	xfi = xati - xfromi
	yfi = yati - yfromi
	zfi = zati - zfromi
c
	fnorm = sqrt(xfi*xfi + yfi*yfi + zfi*zfi)
c
	call normalize(xfi,yfi,zfi)
c
c----- Get up vector
c
	print*,' -- Enter up vector'
	read*, xui, yui, zui
	call normalize(xui,yui,zui)
c
c----- Ensure that initial facing and up vectors are perpendicular
c	Adjust up vector as needed:
c	      u = (f x u) x f
c
	call crossprod(xfi,yfi,zfi,xui,yui,zui,xtemp,ytemp,ztemp)
c
	call crossprod(xtemp,ytemp,ztemp,xfi,yfi,zfi,xui,yui,zui)
c
	call normalize(xui,yui,zui)
c
	print*,'   Inital up vector used - ',xui,yui,zui
	print*,' '
c
	print*,' Enter output choreography filename'
	read(*,'(a)') outfname
c
	open(11,file=outfname,status='unknown')
	rewind(11)
c
c----- Option to output initial frame
c
	print*,' -- Do you want to output the initial frame? (y/n)'
	read(*,'(a)') ans
c
	call lows(ans)
	if (ans .eq. 'y') then
	  write(6,1010) iframe,xfromi,yfromi,zfromi,xati,yati,zati,xui,yui,zui
	  write(11,1010) iframe,xfromi,yfromi,zfromi,xati,yati,zati,xui,yui,zui
 1010	  format(1x,i7,9(1x,1pe11.4))
	endif
c
c=====================================================================
c
c----- Copy initial information to current information
c
	xfromc = xfromi
	yfromc = yfromi
	zfromc = zfromi
c
	xatc = xati
	yatc = yati
	zatc = zati
c
	xuc = xui
	yuc = yui
	zuc = zui
c
c=====================================================================
c
c----- GET INFORMATION ABOUT SHOT
c
  100	print*,' '
	print*,' Enter information about shot to choregraph'
	print*,' '
c
c----- Get number of frames in shot
c
	print*,' -- How many frames in this shot?'
	read*,nframes
c
c----- Choose camera maneuvers for shot
c
	trans_from = .false.
	trans_at = .false.
	circle = .false.
	pan = .false.
	tilt = .false.
	twist = .false.
c
  105	print*,' -- Enter number for camera maneuver'
	print*,'     (0) done selecting maneuvers'
	print*,'     (1) translation of from point'
	print*,'     (2) translation of at point'
	print*,'     (3) circle'
	print*,'     (4) pan'
	print*,'     (5) tilt'
	print*,'     (6) twist'
c
	read*, ichoice
c
	if (ichoice .eq. 0) then
	  goto 500
c
c----- Get parameters for a translation of the from point
c
	elseif (ichoice .eq. 1) then
	  trans_from = .true.
c
  110	  print*,'   -- Which vector to use for from point translation?'
	  print*,'      (1) facing vector at start of shot'
	  print*,'      (2) facing vector at start of frame'
	  print*,'      (3) specify a vector'
	  print*,'      (4) specify a target point & calculate at start of shot'
	  print*,'      (5) specify a target point & calculate at start of frame'
	  read*,num
c
	  trans_from_fc = .false.
	  trans_from_pc = .false.
	  if (num .eq. 1) then
	    trans_from_vec_x = xati - xfromi
	    trans_from_vec_y = yati - yfromi
	    trans_from_vec_z = zati - zfromi
c
	  elseif (num .eq. 2) then
	    trans_from_fc = .true.
c
	  elseif (num .eq. 3) then
	    print*,'     -- Enter vector for from point translation (x y z)'
	    read*, trans_from_vec_x, trans_from_vec_y, trans_from_vec_z
c
	  elseif (num .eq. 4) then
	    print*,'     -- Enter target point for from point translation (x y z)'
	    read*, trans_from_pt_x, trans_from_pt_y, trans_from_pt_z
	    trans_from_vec_x = trans_from_pt_x - xfromi
	    trans_from_vec_y = trans_from_pt_y - yfromi
	    trans_from_vec_z = trans_from_pt_z - zfromi
c
	  elseif (num .eq. 5) then
	    trans_from_pc = .true.
	    print*,'     -- Enter target point for from point translation (x y z)'
	    read*, trans_from_pt_x, trans_from_pt_y, trans_from_pt_z
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 110
	  endif
c
  120	  print*,'   -- How far to translate from point?'
	  print*,'      (1) specify multiple of facing vector at start of shot'
	  print*,'      (2) specify a distance'
	  print*,'      (3) specify a target point and multiple of distance'
	  read*,num
c
	  if (num .eq. 1) then
	    print*,'     Enter multiple of facing vector'
c
	    xtemp = xati - xfromi
	    ytemp = yati - yfromi
	    ztemp = zati - zfromi
	    fnorm = sqrt(xtemp*xtemp + ytemp*ytemp + ztemp*ztemp)
c
	    print*,'       facing vector length = ',fnorm
c
	    read*,zmult
c
	    trans_from_dist_per_frame = zmult*fnorm/float(nframes)
c
	  elseif (num .eq. 2) then
	    print*,'     Enter distance'
	    read*,trans_from_dist
	    trans_from_dist_per_frame = trans_from_dist/float(nframes)
c
	  elseif (num .eq. 3) then
	    print*,'     Enter target point (x y z)'
	    read*, xtemp, ytemp, ztemp
	    print*,'     Enter multiple of distance to target'
	    read*,zmult
c
	    xtemp = xtemp - xfromi
	    ytemp = ytemp - yfromi
	    ztemp = ztemp - zfromi
	    fnorm = sqrt(xtemp*xtemp + ytemp*ytemp + ztemp*ztemp)
c
	    trans_from_dist_per_frame = zmult*fnorm/float(nframes)
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 120
	  endif
c
	  if (.not. trans_from_fc .and. .not. trans_from_pc) then
	    call normalize(trans_from_vec_x,trans_from_vec_y,trans_from_vec_z)
	    trans_from_dx = trans_from_dist_per_frame*trans_from_vec_x
	    trans_from_dy = trans_from_dist_per_frame*trans_from_vec_y
	    trans_from_dz = trans_from_dist_per_frame*trans_from_vec_z
	  endif
c
c----- Get parameters for a translation of the at point
c
	elseif (ichoice .eq. 2) then
	  trans_at = .true.
c
  130	  print*,'   -- Which vector to use for at point translation?'
	  print*,'      (1) facing vector at start of shot'
	  print*,'      (2) facing vector at start of frame'
	  print*,'      (3) specify a vector'
	  print*,'      (4) specify a target point'
	  read*,num
c
	  trans_at_fc = .false.
	  if (num .eq. 1) then
	    trans_at_vec_x = xati - xfromi
	    trans_at_vec_y = yati - yfromi
	    trans_at_vec_z = zati - zfromi
c
	  elseif (num .eq. 2) then
	    trans_at_fc = .true.
c
	  elseif (num .eq. 3) then
	    print*,'     -- Enter vector for at point translation (x y z)'
	    read*, trans_at_vec_x, trans_at_vec_y, trans_at_vec_z
c
	  elseif (num .eq. 4) then
	    print*,'     -- Enter target point for at point translation (x y z)'
	    read*, xtemp, ytemp, ztemp
	    trans_at_vec_x = xtemp - xati
	    trans_at_vec_y = ytemp - yati
	    trans_at_vec_z = ztemp - zati
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 130
	  endif
c
  140	  print*,'   -- How far to translate at point?'
	  print*,'      (1) specify multiple of facing vector at start of shot'
	  print*,'      (2) specify a distance'
	  print*,'      (3) specify a target point and multiple of distance'
	  read*,num
c
	  if (num .eq. 1) then
	    print*,'     Enter multiple of facing vector'
c
	    xtemp = xati - xfromi
	    ytemp = yati - yfromi
	    ztemp = zati - zfromi
	    fnorm = sqrt(xtemp*xtemp + ytemp*ytemp + ztemp*ztemp)
c
	    print*,'        facing vector length = ',fnorm
c
	    read*,zmult
c
	    trans_at_dist_per_frame = zmult*fnorm/float(nframes)
c
	  elseif (num .eq. 2) then
	    print*,'     Enter distance'
	    read*,trans_at_dist
	    trans_at_dist_per_frame = trans_at_dist/float(nframes)
c
	  elseif (num .eq. 3) then
	    print*,'     Enter target point (x y z)'
	    read*, xtemp, ytemp, ztemp
	    print*,'     Enter multiple of distance to target'
	    read*,zmult
c
	    xtemp = xtemp - xfromi
	    ytemp = ytemp - yfromi
	    ztemp = ztemp - zfromi
	    fnorm = sqrt(xtemp*xtemp + ytemp*ytemp + ztemp*ztemp)
c
	    trans_at_dist_per_frame = zmult*fnorm/float(nframes)
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 140
	  endif
c
	  if (.not. trans_at_fc) then
	    call normalize(trans_at_vec_x,trans_at_vec_y,trans_at_vec_z)
	    trans_at_dx = trans_at_dist_per_frame*trans_at_vec_x
	    trans_at_dy = trans_at_dist_per_frame*trans_at_vec_y
	    trans_at_dz = trans_at_dist_per_frame*trans_at_vec_z
	  endif
c
c----- Get parameters for a circle move
c
	elseif (ichoice .eq. 3) then
	  circle = .true.
c
	  print*,'   -- Enter degrees to circle for shot'
	  print*,'      (counterclockwise is positive)'
	  read*,circle_degrees
	  circle_degrees_per_frame = circle_degrees/float(nframes)
c
  150	  print*,'   -- Which point to circle around?'
	  print*,'      (1) at point at start of shot'
	  print*,'      (2) at point at start of frame'
	  print*,'      (3) specify a point'
	  read*,num
c
	  circle_atc = .false.
	  if (num .eq. 1) then
	    circle_pt_x = xati
	    circle_pt_y = yati
	    circle_pt_z = zati
c
	  elseif (num .eq. 2) then
	    circle_atc = .true.
c
	  elseif (num .eq. 3) then
	    print*,'     -- Enter point to circle (x y z)'
	    read*, circle_pt_x, circle_pt_y, circle_pt_z
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 150
	  endif
c
  160	  print*,'   -- Which vector to use as circle rotation axis?'
	  print*,'      (1) up vector at start of shot'
	  print*,'      (2) up vector at start of frame'
	  print*,'      (3) specify a vector'
	  read*,num
c
	  circle_uc = .false.
	  if (num .eq. 1) then
	    circle_vec_x = xui
	    circle_vec_y = yui
	    circle_vec_z = zui
c
	  elseif (num .eq. 2) then
	    circle_uc = .true.
c
	  elseif (num .eq. 3) then
	    print*,'     -- Enter vector for circle rotation axis (x y z)'
	    read*, circle_vec_x, circle_vec_y, circle_vec_z
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 160
	  endif
c
c----- Get parameters for a camera pan move
c
	elseif (ichoice .eq. 4) then
	  pan = .true.
c
	  print*,'   -- Enter degrees to pan for shot'
	  print*,'      (counterclockwise is positive)'
	  read*,pan_degrees
	  pan_degrees_per_frame = pan_degrees/float(nframes)
c
  210	  print*,'   -- Which vector to use as pan rotation axis?'
	  print*,'      (1) up vector at start of shot'
	  print*,'      (2) up vector at start of frame'
	  read*,num
c
	  pan_uc = .false.
	  if (num .eq. 1) then
	    pan_vec_x = xui
	    pan_vec_y = yui
	    pan_vec_z = zui
c
	  elseif (num .eq. 2) then
	    pan_uc = .true.
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 210
	  endif
c
c----- Get parameters for a camera tilt move
c
	elseif (ichoice .eq. 5) then
	  tilt = .true.
c
	  print*,'   -- Enter degrees to tilt for shot'
	  print*,'      (counterclockwise is positive)'
	  read*,tilt_degrees
	  tilt_degrees_per_frame = tilt_degrees/float(nframes)
c
  260	  print*,'   -- Which vector to use as tilt rotation axis?'
	  print*,'      (1) side vector at start of shot'
	  print*,'      (2) side vector at start of frame'
	  read*,num
c
	  tilt_sc = .false.
	  if (num .eq. 1) then
	    xfi = xati - xfromi
	    yfi = yati - yfromi
	    zfi = zati - zfromi
	    call crossprod(xfi,yfi,zfi,xui,yui,zui,
     &	                   tilt_vec_x,tilt_vec_y,tilt_vec_z)
c
	  elseif (num .eq. 2) then
	    tilt_uc = .true.
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 260
	  endif
c
c----- Get parameters for a camera twist move
c
	elseif (ichoice .eq. 6) then
	  twist = .true.
c
	  print*,'   -- Enter degrees to twist for shot'
	  print*,'      (counterclockwise is positive)'
	  read*,twist_degrees
	  twist_degrees_per_frame = twist_degrees/float(nframes)
c
  310	  print*,'   -- Which vector to use as twist rotation axis?'
	  print*,'      (1) facing vector at start of shot (negated)'
	  print*,'      (2) facing vector at start of frame (negated)'
	  read*,num
c
	  twist_fc = .false.
	  if (num .eq. 1) then
	    twist_vec_x = xfromi - xati
	    twist_vec_y = yfromi - yati
	    twist_vec_z = zfromi - zati
c
	  elseif (num .eq. 2) then
	    twist_fc = .true.
c
	  else
	    print*,'     BAD INPUT'
	    print*,' '
	    goto 310
	  endif
c
c----- Handle bad choice of camera maneuver
c
	else
	  print*,'     BAD INPUT'
	  print*,' '
	endif
c
c----- Loop back to camera maneuver choice
c
	goto 105
c
c=====================================================================
c
c----- Copy current information to frame information
c
  500	xfromf = xfromc
	yfromf = yfromc
	zfromf = zfromc
c
	xatf = xatc
	yatf = yatc
	zatf = zatc
c
	xff = xfc
	yff = yfc
	zff = zfc
c
	xuf = xuc
	yuf = yuc
	zuf = zuc
c
c----- Loop over frames in shot
c
	do jf = 1, nframes
c
c----- Choreograph translation of from point
c
	  if (trans_from) then
c
	    if (trans_from_fc) then
	      trans_from_vec_x = xatc - xfromc
	      trans_from_vec_y = yatc - yfromc
	      trans_from_vec_z = zatc - zfromc
c
	      call normalize(trans_from_vec_x,trans_from_vec_y,trans_from_vec_z)
	      trans_from_dx = trans_from_dist_per_frame*trans_from_vec_x
	      trans_from_dy = trans_from_dist_per_frame*trans_from_vec_y
	      trans_from_dz = trans_from_dist_per_frame*trans_from_vec_z
c
	    elseif (trans_from_pc) then
	      trans_from_vec_x = trans_from_pt_x - xfromc
	      trans_from_vec_y = trans_from_pt_y - yfromc
	      trans_from_vec_z = trans_from_pt_z - zfromc
c
	      call normalize(trans_from_vec_x,trans_from_vec_y,trans_from_vec_z)
	      trans_from_dx = trans_from_dist_per_frame*trans_from_vec_x
	      trans_from_dy = trans_from_dist_per_frame*trans_from_vec_y
	      trans_from_dz = trans_from_dist_per_frame*trans_from_vec_z
	    endif
c
	    xfromf = xfromf + trans_from_dx
	    yfromf = yfromf + trans_from_dy
	    zfromf = zfromf + trans_from_dz
	  endif
c
c----- Choreograph translation of at point
c
	  if (trans_at) then
c
	    if (trans_at_fc) then
	      trans_at_vec_x = xatc - xfromc
	      trans_at_vec_y = yatc - yfromc
	      trans_at_vec_z = zatc - zfromc
c
	      call normalize(trans_at_vec_x,trans_at_vec_y,trans_at_vec_z)
	      trans_at_dx = trans_at_dist_per_frame*trans_at_vec_x
	      trans_at_dy = trans_at_dist_per_frame*trans_at_vec_y
	      trans_at_dz = trans_at_dist_per_frame*trans_at_vec_z
	    endif
c
	    xatf = xatf + trans_at_dx
	    yatf = yatf + trans_at_dy
	    zatf = zatf + trans_at_dz
	  endif
c
c----- Choreograph circle move
c
	  if (circle) then
c
	    if (circle_atc) then
	      circle_pt_x = xatc
	      circle_pt_y = yatc
	      circle_pt_z = zatc
	    endif
c
	    if (circle_uc) then
	      circle_vec_x = xuc
	      circle_vec_y = yuc
	      circle_vec_z = zuc
	    endif
c
c----- Rotate from point
c
	    x1 = circle_pt_x - xfromf
	    y1 = circle_pt_y - yfromf
	    z1 = circle_pt_z - zfromf
c
	    call rotate(x1,y1,z1,circle_vec_x,circle_vec_y,circle_vec_z,
     &	                circle_degrees_per_frame,x2,y2,z2)
c
	    xfromf = circle_pt_x - x2
	    yfromf = circle_pt_y - y2
	    zfromf = circle_pt_z - z2
c
c----- Rotate up vector
c
	    x1 = xuf
	    y1 = yuf
	    z1 = zuf
c
	    call rotate(x1,y1,z1,circle_vec_x,circle_vec_y,circle_vec_z,
     &	                circle_degrees_per_frame,x2,y2,z2)
c
	    xuf = x2
	    yuf = y2
	    zuf = z2
c
	  endif
c
c----- Choreograph camera pan move
c
	  if (pan) then
c
	    if (pan_uc) then
	      pan_vec_x = xuc
	      pan_vec_y = yuc
	      pan_vec_z = zuc
	    endif
c
	    x1 = xatf - xfromf
	    y1 = yatf - yfromf
	    z1 = zatf - zfromf
c
	    call rotate(x1,y1,z1,pan_vec_x,pan_vec_y,pan_vec_z,
     &	                pan_degrees_per_frame,x2,y2,z2)
c
	    xatf = x2 + xfromf
	    yatf = y2 + yfromf
	    zatf = z2 + zfromf
c
	  endif
c
c----- Choreograph camera tilt move
c
	  if (tilt) then
c
	    if (tilt_uc) then
	      x1 = xatc - xfromc
	      y1 = yatc - yfromc
	      z1 = zatc - zfromc
c
	      call crossprod(x1,y1,z1,xuc,yuc,zuc,
     &	                     tilt_vec_x,tilt_vec_y,tilt_vec_z)
	    endif
c
	    x1 = xatf - xfromf
	    y1 = yatf - yfromf
	    z1 = zatf - zfromf
c
	    call rotate(x1,y1,z1,tilt_vec_x,tilt_vec_y,tilt_vec_z,
     &	                tilt_degrees_per_frame,x2,y2,z2)
c
	    xatf = x2 + xfromf
	    yatf = y2 + yfromf
	    zatf = z2 + zfromf
c
	  endif
c
c----- Choreograph camera twist move
c
	  if (twist) then
c
	    if (twist_fc) then
	      twist_vec_x = xfromc - xatc
	      twist_vec_y = yfromc - yatc
	      twist_vec_z = zfromc - zatc
	    endif
c
	    call rotate(xuf,yuf,zuf,twist_vec_x,twist_vec_y,twist_vec_z,
     &	                twist_degrees_per_frame,x2,y2,z2)
c
	    xuf = x2
	    yuf = y2
	    zuf = z2
c
	  endif
c
c----- Write out information for this frame
c
	  iframe = iframe + 1
	  write(6,1010) iframe,xfromf,yfromf,zfromf,xatf,yatf,zatf,xuf,yuf,zuf
	  write(11,1010) iframe,xfromf,yfromf,zfromf,xatf,yatf,zatf,xuf,yuf,zuf
c
c----- Copy frame information to current information
c
	  xfromc = xfromf
	  yfromc = yfromf
	  zfromc = zfromf
c
	  xatc = xatf
	  yatc = yatf
	  zatc = zatf
c
	  xfc = xff
	  yfc = yff
	  zfc = zff
c
	  xuc = xuf
	  yuc = yuf
	  zuc = zuf
c
	enddo
c
c=====================================================================
c
	print*,' '
	print*,' Another shot?'
	read(*,'(a)') ans
c
	call lows(ans)
	if (ans .eq. 'y') then
c
c----- Copy current information to initial information
c
	  xfromi = xfromc
	  yfromi = yfromc
	  zfromi = zfromc
c
	  xati = xatc
	  yati = yatc
	  zati = zatc
c
	  xfi = xfc
	  yfi = yfc
	  zfi = zfc
c
	  xui = xuc
	  yui = yuc
	  zui = zuc
c
	  goto 100
c
	endif
c
  999	close(11)
c
	end