;
;	Procedure InSARFringe_wrapped.PRO
;
;	This procedure plots a map of the faults for the model
;	seismicity generated by VC_SER.f, then computes and plots
;	the equivalent, time-dependent wrapped InSAR Fringes
;	assuming either C or L band (user chooses).
;	
;	We sometimes use the Massonet et al. look vector as standard:
;	(x,y,z) = (.33,.07,.94)
;
;	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 insarfringe_wrapped, 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)

	InSAR_band = 'C'

	
	
	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

	NX = long(sqrt(npoint))
	NY = long(sqrt(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
	
;	m = long(-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
;	m = m+1
	xpt(i) = xpt_read(i)
	ypt(i) = ypt_read(i)
	dux(i) = dux_read(i)
	duy(i) = duy_read(i)
;	endelse
	endfor
	
;	npoint = m
	
	xpt_center_1	=	fltarr(npoint,/no)
	xpt_center_2	=	fltarr(npoint,/no)
	ypt_center_1	=	fltarr(npoint,/no)
	ypt_center_2	=	fltarr(npoint,/no)
	
	fring_val_cum 	=	fltarr(npoint,/no)
	fring_val_int 	=	fltarr(npoint,/no)
	duz 		=	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
	
;	fringe_val_cum is the dot product of the look vector with the
;		deformation vector

	print, ' '
	print, 'Enter C- or L-band frequency for InSAR:  (C,L)'
	read, InSAR_band

;
;	Note that L-band is at 1.3 GHz

	lambda = 5.7
	if (InSAR_band eq 'L') then lambda = 23.07


;	fringe_val_cum = (.33*dux + .07*duy + .94*duz)* vfact/lambda
;
;	Change look angle to horizontal, from SE
;

	look_x = .707
	look_y = .707
	look_z = 0.

	fringe_val_cum = (look_x*dux + look_y*duy + look_z*duz)* vfact/lambda



	X_ARRAY		=	fltarr(NX,NY,/no)
	Y_ARRAY		= 	fltarr(NY,NY,/no)
	FRINGE_VAL	=	fltarr(NX,NY,/no)


	m = -1
	for i=0L,NX-1 do begin
		for j = 0L, NY-1 do begin
		m = m + 1
		X_ARRAY(i,j) = xpt_center_1(m)
		Y_ARRAY(i,j) = ypt_center_1(m)
		FRINGE_VAL(i,j) = fringe_val_cum(m)
		endfor
	endfor

	fringe_int = 1.


	fringe_inv = 1./fringe_int

	fringe_cum_max = max(fringe_val_cum)
	fringe_cum_min = min(fringe_val_cum)

	print, ' Max, Min fringe values:  ', fringe_cum_max, fringe_cum_min

	fringe_cum_diff = fringe_cum_max - fringe_cum_min

	fringe_min = float(fix(fringe_cum_min * fringe_inv))/fringe_inv


	fringe_n = (fringe_cum_max - fringe_cum_min)*fringe_inv

	num_fringe = fix(abs(fringe_n + 1))
	if(num_fringe gt 56) then num_fringe = 56

;
;	Find the maximum and minimum for the cumulative fringe
;		values


	
	xpt_center_2 = xpt_center_1 + dux
	ypt_center_2 = ypt_center_1 + duy
	
	

;

      	 	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

	rebin_factor = 15
	
	NEW_DIM_X = rebin_factor * NX
	NEW_DIM_Y = rebin_factor * NY

	NEW_FRINGE_VAL  =  REBIN(FRINGE_VAL, NEW_DIM_X, NEW_DIM_Y)
	NEW_X_ARRAY  = REBIN(X_ARRAY, NEW_DIM_X, NEW_DIM_Y)
	NEW_Y_ARRAY  = REBIN(Y_ARRAY, NEW_DIM_X, NEW_DIM_Y)

	TWO_PI  = 6.283185307

	NEW_FRINGE_VAL =  sin (TWO_PI * NEW_FRINGE_VAL)

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

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

	if (InSAR_band eq 'C') then begin

     plot, [0],[0], $
       xstyle=1,ystyle=1, $              ; extend axes
	xticks = 1, yticks = 1, 	$
        xrange=[xleft,xright], $
        yrange=[ydown,yup], $
        title='WRAPPED SYNTHETIC InSAR FRINGES:  C - BAND'

	endif else begin

     plot, [0],[0], $
       xstyle=1,ystyle=1, $              ; extend axes
	xticks = 1, yticks = 1, 	$
        xrange=[xleft,xright], $
        yrange=[ydown,yup], $
        title='WRAPPED SYNTHETIC InSAR FRINGES:  L - BAND'


	endelse

	loadct, 13	; rainbow


	contour, NEW_FRINGE_VAL, NEW_X_ARRAY, NEW_Y_ARRAY,  $
	levels = findgen(60) * .03333333333 - 1.,   $
	/FILL, /OVERPLOT, c_colors = findgen(60)*220/60 + 35

	line_thick = 6.
	doub_line_thick = 18.


;
;	Add descriptive words
;

	xs0 = -400
	ys0 = -400.
	xs1 = xs0
	ys1 = -430.
	xs2 = 0.
	ys2 = -400.
;	ys2 = - .65 * AXLEN


;	astr0 = 'Average Velocity Over Years:'
	astr0 = 'Displacement Over Years:'
	astr1 = strcompress(string(tbase)) + ' - ' + strcompress(string(tlater))
	astr2 = 'Contour Interval = '
	astr3 = ' * Wavelength'
;	astr2 = astr2 + strcompress(string(fringe_int)) + astr3
	astr2 = astr2 + strmid(strcompress(string(fringe_int)),0,4) + astr3
	xyouts, xs0, ys0, astr0, charsize = 1.00
	xyouts, xs1, ys1, astr1, charsize = 1.00
	xyouts, xs2, ys2, astr2, charsize = 1.00	

;
;	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 = line_thick, color = 0



	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 = line_thick, color = 0

	x1(0) = XR(nmaxx)
	y1(0) = YR(nmaxx)
	x1(1) = XBR
	y1(1) = YBR
	oplot, x1, y1, linestyle = 2, thick = line_thick, color = 0
;
;	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=doub_line_thick, color = 0
	endif
	endfor


	astr3  =  'Look'
	astr4  =  'Vector: '

	xs3 = 390.
	ys3 = 375.
	xs4 = 390.
	ys4 = 350.
	xyouts, xs3, ys3, astr3, charsize = 1.0
	xyouts, xs4, ys4, astr4, charsize = 1.0

	x_look_left = 400.
	y_look_left = 300.
	x_look_right = x_look_left + 30. * look_x
	y_look_right = y_look_left + 30. * look_y
	arrow,x_look_left, y_look_left, x_look_right, y_look_right, $
		thick = 2., /data


	


;
	
;

;
;
       end