;
;	Procedure MODEL.PRO
;
;	This procedure plots a 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
;	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, nfault, vplatx, vplaty, xfe, xfw, yfe, yfw
;
;	****************************************************************
;
	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)


     	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
;
;	Find whether plate boundary goes out top or bottom
;
;    	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

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


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

;
;	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)


;	oplot, x1, y1, psym = 8, symsize = .7, thick = 8.0, color = 100
	oplot, x1, y1, linestyle = 0, thick = 8.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 = 8.0, color = 100

	x1(0) = XR(nmaxx)
	y1(0) = YR(nmaxx)
	x1(1) = XBR
	y1(1) = YBR
	oplot, x1, y1, linestyle = 2,thick = 8.0, color = 100
;
;
       end