;
;	Procedure DISPL_ARROW.PRO
;
;	Plots a map of horizontal deformation arrows superposed on
;		a fault map.
;
;	Primary Variables:
;
;
;	tstep 		- Size of time steps in yr.
;	nfault 		- number of fault segments
;	vplatx		- East-component of plate velocity in cm/yr
;	vplaty		- North-component of plate velocity in cm/yr
;	xfe(i) 		- x-component of easternmost point of fault segment i
;	yfe(i) 		- y-component of easternmost point of fault segment i
;	xfw(i) 		- x-component of westernmost point of fault segment i
;	yfw(i) 		- y-component of westernmost point of fault segment i
;	tslip(i,j) 	- Time at which segment i slipped for jth time (yr)
;	islip(i)   	- Number of times segment i slipped (# of earthquakes)

;
;	****************************************************************
;
	pro displ_arrow, tstep, nfault, vplatx, vplaty,$
		xfe, xfw, yfe, yfw, islip, tslip
;
;	****************************************************************
;
	NF = long(nfault)
;
	XL	=	fltarr(NF,/no)
	XR	=	fltarr(NF,/no)
	YL	=	fltarr(NF,/no)
	YR	=	fltarr(NF,/no)
	XINT	=	fltarr(NF,/no)
	YINT	=	fltarr(NF,/no)
	
	
	fname = 'deform_data.d'
;	print, 'Enter Name of File with DISPLACEMENT VECTORS: '
;	read, fname
;	
	openr, 20, fname
	readf, 20, vplate, vfact, tbase, tlater, tmin, npoint
	
	xpt		=	fltarr(npoint,/no)
	ypt		=	fltarr(npoint,/no)
	dux		=	fltarr(npoint,/no)
	duy		=	fltarr(npoint,/no)
	xpt_read	=	fltarr(npoint,/no)
	ypt_read	=	fltarr(npoint,/no)
	dux_read	=	fltarr(npoint,/no)
	duy_read	=	fltarr(npoint,/no)
	
	readf, 20, xpt_read
	readf, 20, ypt_read
	readf, 20, dux_read
	readf, 20, duy_read
	
	close, 20
	


;
;	Locate east, west, north, south limits on faults
;
	xmin1 = min(xfe, nminx1)
	xmin2 = min(xfw, nminx2)
	xmin = xmin1
	nminx = nminx1
	if (xmin2 lt xmin1) then begin
	xmin = xmin2
	nminx = nminx2
	endif

;
	xmax1 = max(xfe, nmaxx1)
	xmax2 = max(xfw, nmaxx2)
	xmax = xmax1
	nmaxx = nmaxx1
	if (xmax2 gt xmax1) then begin
	xmax = xmax2
	nmaxx = nmaxx2
	endif 

;
	ymin1 = min(yfe, nminy1)
	ymin2 = min(yfw, nminy2)
	ymin = ymin1
	nminy = nminy1

	if (ymin2 lt ymin1) then begin
	ymin = ymin2
	nminy = nminy2

	endif

;
	ymax1 = max(yfe, nmaxy1)
	ymax2 = max(yfw, nmaxy2)
	ymax = ymax1
	nmaxy = nmaxy1

	if (ymax2 gt ymax1) then begin
	ymax = ymax2
	nmaxy = nmaxy2

	endif
;

	resp2 = 'y'

;
;	Find Center of Picture
;
	xcenter = (xmin + xmax)/2.
	ycenter = (ymin + ymax)/2.
;
	xlong = 1.25*(xcenter - xmin)
	ylong = 1.25*(ycenter - ymin)


     	axlen = max(xlong,ylong)
	xleft  = -axlen
	xright =  axlen
	
	yup   =  axlen
	ydown = -axlen
	
	xpt_read = xpt_read - xcenter
	ypt_read = ypt_read - ycenter
	
	l = -1
	for i = 0L,npoint-1 do begin
	if (xpt_read(i) lt 0. and $
		ypt_read(i) lt -.6*axlen) then begin
;
;	We don't want these points
;
	endif else begin
	l = l+1
	xpt(l) = xpt_read(i)
	ypt(l) = ypt_read(i)
	dux(l) = dux_read(i) * vfact
	duy(l) = duy_read(i) * vfact
	endelse
	endfor
	
	npoint = l
	
	xpt_center_1	=	fltarr(npoint,/no)
	xpt_center_2	=	fltarr(npoint,/no)
	ypt_center_1	=	fltarr(npoint,/no)
	ypt_center_2	=	fltarr(npoint,/no)

	du_max 		=	fltarr(npoint,/no)
	
;	xpt_center_1 = xpt - xcenter
;	ypt_center_1 = ypt - ycenter

	for i = 0L, npoint-1 do begin
	xpt_center_1(i) = xpt(i)
	ypt_center_1(i) = ypt(i)
	endfor


	du_scale = 50.

	resp_S = 'V'
	print, ' '
	print, ' Scale to Plate Velocity (V) or '
	print, ' Max Displacement (D) ?'
	read, resp_S

	if (resp_S eq 'D') then begin

	for i = 0L,npoint-1 do begin
	du_max(i) = sqrt(dux(i)*dux(i) + duy(i)*duy(i)); put CM units back in
	endfor

	max_displ = max(du_max)
	

	print, 'max_displ=', max_displ

;	Test du_max to see how big they are.  Keep only those that are
;	large enough

	dux(*) = dux(*) / max_displ
	duy(*) = duy(*) / max_displ

	dux(*) =  dux(*) * 5. * du_scale
	duy(*) =  duy(*) * 5. * du_scale

	xpt_center_2 = xpt_center_1 + dux
	ypt_center_2 = ypt_center_1 + duy

	endif else begin

	dux(*) =  dux(*) * (du_scale * .75)/(.5*vfact)
	duy(*) =  duy(*) * (du_scale * .75)/(.5*vfact)



	xpt_center_2 = xpt_center_1 + dux
	ypt_center_2 = ypt_center_1 + duy

	endelse
	

	
;

      	 	XL = xfw - xcenter
    		XR = xfe - xcenter
       		YL = yfw - ycenter
       		YR = yfe - ycenter
      		XINT = .5*(XL + XR)
      		YINT = .5*(YL + YR)

;
      slope= -vplaty/vplatx
;
;	Find whether plate boundary goes out top or bottom
;
;    	Do left boundary first
;


	B = YL(nminx) - slope*XL(nminx)
	B1 = slope*XL(nminx) + B


	if (B1 gt AXLEN) then begin
;    It goes out the top
	YBL = AXLEN
	XBL = (YBL - B1)/slope
	endif else begin
;    It goes out the side
	XBL= -AXLEN
	YBL= slope*XBL + B
	endelse

	XBL = 10.* XBL
	YBL = slope*XBL + B

;    Now do right boundary

	B = YR(nmaxx) - slope*XR(nmaxx)
	C = slope*AXLEN + B

	if (C lt -AXLEN) then begin
;    It goes out the bottom
	YBR = -AXLEN
	XBR = -B/slope
	endif else begin
;    It goes out the side
	XBR = AXLEN
	YBR = C
	endelse

	XBR = 10.*XBR
	YBR = slope*XBR + B

   iflag = intarr(nfault)

   for j = 0L, nfault-1 do begin
     iflag(j) = n_elements(where(sqrt((XINT(j)-XINT)^2 $
                                    + (YINT(j)-YINT)^2) lt 1.e-5))-1
   endfor

	!P.POSITION =  [0., 0., 1., 1.]

	xleft = -450.
	xright = 450.
	ydown = -450.
	yup = 450.

	loadct, 13	; rainbow

   plot,[0],[0], $                          ; dummies
        xstyle=1,ystyle=1, $              ; extend axes
        xrange=[xleft,xright], $           
        yrange=[ydown,yup], $
	xticks= 1, yticks = 1, $
        title='Model Faults:  Map View', color = 0

;	
;	Add descriptive words
;

        
	xs0 = -1.15*AXLEN
	ys0 = -.75*AXLEN
	xs1 = xs0
	ys1 = -.85*AXLEN
	xs2 = xs0
	ys2 = - 1.0 * AXLEN


	if(resp_S eq 'V') then begin
	astr0 = 'Average Velocity Over Years:'
	astr1 = strcompress(string(tbase)) + $
		' - ' + strcompress(string(tlater))
	astr2 = '17.5 mm/yr: '
	xyouts, xs0, ys0, astr0, charsize = 1.0
	xyouts, xs1, ys1, astr1, charsize = 1.0
	xyouts, xs2, ys2, astr2, charsize = 1.0

	endif else begin

	scale_du  =  max_displ * .2
	astr0 = 'Total Displacement Over Years:'
	astr1 = strcompress(string(tbase)) $
		+ ' - ' + strcompress(string(tlater))
	astr2 = 'Scale = ' + strmid(strcompress(string(scale_du)),0,4) + ' CM: '
	xyouts, xs0, ys0, astr0, charsize = 1.0
	xyouts, xs1, ys1, astr1, charsize = 1.0
	xyouts, xs2, ys2, astr2, charsize = 1.0

	endelse	

;
;	Loop over faults
;
;      !P.charsize = 1.5

;	A = findgen(16) * (!pi*2/16.)

;	usersym, cos(A), sin(A), /fill

	x1	=	fltarr(2,/no)
	y1	=	fltarr(2,/no)

	for i=0L,NF-1 do begin
	x1(0) = XL(i)
	y1(0) = YL(i)
	x1(1) = XR(i)
	y1(1) = YR(i)
;
;	Do not plot end points on fault segments
;
;	oplot, x1, y1, psym = 8, symsize = .7
	oplot, x1, y1, linestyle = 0, thick = 6.0, color = 100

	endfor
;
;	Plot plate motion extenders
;

	x1(0) = XL(nminx)
	y1(0) = YL(nminx)
	x1(1) = XBL
	y1(1) = YBL
	oplot, x1, y1, linestyle = 2, thick = 6.0, color = 100

	x1(0) = XR(nmaxx)
	y1(0) = YR(nmaxx)
	x1(1) = XBR
	y1(1) = YBR
	oplot, x1, y1, linestyle = 2, thick = 6.0, color = 100
;
;	Draw faults that slipped

	time_mat = fltarr(NF,/no)

  	for fault_index = 0L,NF-1 do begin
	time_mat(fault_index) = 0
	endfor
	
 	for fault_index = 0L,NF-1 do begin
		for time_index = 0L,islip(fault_index) -1 do begin
		time_test_1 = tslip(fault_index,time_index)- tstep*.001
		time_test_2 = tslip(fault_index,time_index)- tstep*.001

		if (time_test_1 gt tbase and time_test_2 le tlater) then begin
		time_mat(fault_index) = 1
		endif
		
		endfor
   	endfor
;
;	Draw slipped segments
;
	for i=0L,NF-1 do begin
	if (time_mat(i) eq 1) then begin
	x1(0) = XL(i)
	y1(0) = YL(i)
	x1(1) = XR(i)
	y1(1) = YR(i)
	oplot, x1, y1, linestyle=0, thick=14.0, color = 150
	endif
	endfor


;
;	Plot Arrows
;

	xscale_left = -225.
	yscale_left = ys2 + 10.

	if (resp_S eq 'D') then begin
	xscale_right = xscale_left + (du_scale*5.*.2)
	endif else begin
	xscale_right = xscale_left + (du_scale * .75)
	endelse

	yscale_right = yscale_left

	for i = 0L, npoint-1 do begin
	
	arrow, xpt_center_1, ypt_center_1, xpt_center_2, ypt_center_2,	$
		/data, color = 254


	endfor

	arrow, xscale_left, yscale_left, xscale_right, yscale_right, $
		/data, color = 254
		
		
;

;
;
       end