;
;	Procedure MODEL_3D.PRO
;
;	This procedure plots a three-dimensional map of the faults
;	for the model seismicity generated by VC_SER.f
;
;	Primary Variables:
;
;	nfault 		- number of fault segments
;	vplatx		- East-component of plate velocity in cm/yr
;	vplaty		- North-component of plate velocity in cm/yr
;	db(i)		- Depth to bottom, segment i (km)
;	dt(i)		- Depth to top, segment i (km)
;	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
;
;	****************************************************************
;
	pro model_3d, nfault, vplatx, vplaty, xfe, xfw, yfe, yfw,	$
		db, dt
;
;	****************************************************************
;
	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)

;
;	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
;
;	Print Extremal Patches:
;
	print,' '
	print,' Westernnmost patch is: '
	print, 'Patch:  ', nminx
	print, 'xfe,yfe,xfw,yfw:  ', xfe(nminx), yfe(nminx), xfw(nminx),$
			 		yfw(nminx)
;
	print,' '
	print,' Easternnmost patch is: '
	print, 'Patch:  ', nmaxx
	print, 'xfe,yfe,xfw,yfw:  ', xfe(nmaxx), yfe(nmaxx), xfw(nmaxx),$
					yfw(nmaxx)

	print,' '
	print,' Northernnmost patch is: '
	print, 'Patch:  ', nmaxy
	print, 'xfe,yfe,xfw,yfw:  ', xfe(nmaxy), yfe(nmaxy), xfw(nmaxy),$
					yfw(nmaxy)

	print,' '
	print,' Southernnmost patch is: '
	print, 'Patch:  ', nminy
	print, 'xfe,yfe,xfw,yfw:  ', xfe(nminy), yfe(nminy), xfw(nminy),$
					yfw(nminy)

;

	resp2 = 'y'
      	print,' '
;      	print,' Omit numbers on the faults? (y/n)'
;      	read, format='(a1)', resp2

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

	xlong = .9*(xcenter - xmin)
	ylong = .9*(ycenter - ymin)



     	axlen = max(xlong,ylong)
	xleft  = -axlen
	xright =  axlen
	
	yup   =  axlen
	ydown = -axlen

;

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

;
      slope= -vplaty/vplatx
;
;	No plate motion extenders on 3-D plot
;
;    	Do left boundary first
;


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

	print, 'B1', B1

	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

	print, 'B,C', B,C

	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

;	loadct, 13	; rainbow
	loadct, 39	; rainbow + white

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



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

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

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

	x1	=	fltarr(5)
	y1	=	fltarr(5)
	z1	=	fltarr(5)

	x2 = fltarr(2)
	y2 = fltarr(2)
	z2 = fltarr(2)

	num_major_faults = long(12)

	for i=0L,NF-1 do begin

	if (i ge 0 and i le 61) then jfault = 1
	if (i ge 62 and i le 90) then jfault = 2
	if (i ge 91 and i le 114) then jfault = 3
	if (i ge 115 and i le 132) then jfault = 4
	if (i ge 133 and i le 142) then jfault = 5
	if (i ge 143 and i le 166) then jfault = 6
	if (i ge 167 and i le 173) then jfault = 7
	if (i ge 174 and i le 178) then jfault = 8
	if (i ge 179 and i le 188) then jfault = 9
	if (i ge 189 and i le 193) then jfault = 10
	if (i ge 194 and i le 204) then jfault = 11
	if (i ge 205 and i le 214) then jfault = 12



;	color_value = fix(float(jfault) *  220./float(num_major_faults)   +  20.)
	color_value = 60

	x1(0) = XL(i)
	y1(0) = YL(i)
	z1(0) = dt(i)

	x1(1) = XR(i)
	y1(1) = YR(i)
	z1(1) = dt(i)

	x1(2) = XR(i)
	y1(2) = YR(i)
	z1(2) = -1.d0*db(i)

	x1(3) = XL(i)
	y1(3) = YL(i)
	z1(3) = -1.d0*db(i)

	x1(4) = XL(i)
	y1(4) = YL(i)
	z1(4) = dt(i)

	x2(0) = XL(i)
	y2(0) = YL(i)
	z2(0) = 0.

	x2(1) = XR(i)
	y2(1) = YR(i)
	z2(1) = 0.


   	plot_3dbox, x1, y1, z1, $                          ; dummies
        xstyle=4,ystyle=4, zstyle = 4,$              ; surpress axes
        xrange=[xleft,xright], $
        yrange=[ydown,yup], $
	zrange = [-40.,20.],	$
	ax = 85., az = 10.,	$
        title='Model Faults:  Map View', /noerase, linestyle = 0, $
;	thick = 4.0, color = 100
	thick = 2., color = color_value

   	plot_3dbox, x2, y2, z2, $                          ; dummies
        xstyle=4,ystyle=4, zstyle = 4,$              ; surpress axes
        xrange=[xleft,xright], $
        yrange=[ydown,yup], $
	zrange = [-40.,20.],	$
	ax = 85., az = 10.,	$
        title='Model Faults:  Map View', /noerase, linestyle = 0, $
;	thick = 4.0, color = 100
	thick = 8., color = 0

	xs = -.975*axlen
	ys =  -.9*axlen
;	astr = 'Color-Coded Fault Model'
	astr = 'Fault Model'
	xyouts, xs, ys, z=0., astr, charsize = 1.5, color = 0	


	endfor

;	if (num_major_faults eq 12) then begin
	if (num_major_faults eq 0) then begin
	for i = 0L, NF-1 do begin

	if (i eq 30) then begin
	xs = XL(i)
	ys = YL(i)- 55.
	astr = 'SA'
	xyouts, xs, ys, z = 0., astr, charsize = 1.15 , color = 0
	endif

	if (i eq 62) then begin
	xs = XL(i)
	ys = YL(i)- 32.
	astr = 'SJ'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 100) then begin
	xs = XL(i)
	ys = YL(i)- 35.
	astr = 'ELS'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 122) then begin
	xs = XL(i)
	ys = YL(i)
	astr = 'IV'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 135) then begin
	xs = XL(i)
	ys = YL(i)
	astr = 'LS'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 160) then begin
	xs = XL(i)
	ys = YL(i)  -65.
	astr = 'GAR'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 170) then begin
	xs = XL(i)
	ys = YL(i)  -50.
	astr = 'PV'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 178) then begin
	xs = XL(i)
	ys = YL(i)  -65.
	astr = 'SCI'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 180) then begin
	xs = XL(i)  + 40.
	ys = YL(i)  -70.
	astr = 'PISG'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 190) then begin
	xs = XL(i)
	ys = YL(i) + 30.
	astr = 'BR'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 196) then begin
	xs = XL(i) - 10.
	ys = YL(i)  -57.5
	astr = 'SM'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif

	if (i eq 212) then begin
	xs = XL(i) - 55.
	ys = YL(i)  -30.
	astr = 'LAN'
	xyouts, xs, ys, z = 0., astr, charsize = 1. , color = 0
	endif
	endfor	
	endif

	x1(0) = 1000.
	y1(0) = 1000.
	z1(0) = 0.

	x1(1) = 1000.
	y1(1) = 1000.
	z1(1) = 0.

	x1(2) = 1000.
	y1(2) = 1000.
	z1(2) = 0.

	x1(3) = 1000.
	y1(3) = 1000.
	z1(3) = 0.

	x1(4) = 1000.
	y1(4) = 1000.
	z1(4) = 0.

;	Get rid of those annoyingly weird axes left over, the stuff I don't want by
;		overplotting them in white!!!

  	plot_3dbox, x1, y1, z1, $                          ; dummies
        xstyle=4,ystyle=4, zstyle = 4,$              ; surpress axes
        xrange=[xleft,xright], $
        yrange=[ydown,yup], $
	zrange = [-40.,20.],	$
	ax = 85., az = 10.,	$
        title='Model Faults:  Map View', /noerase, linestyle = 0, $
	thick = 8.0, color = 255
;
;
       end