FUNCTION mergeline, initinput, ptopass, extras, anglethreshold=anglethreshold ;extras=[Bms,Bearth,num]

;this function takes the starting point and increments in the j direction
;in steps of 0.16Re both sides of the merging line until the the threshold
;test fails or until the maximum imposed length is reached (controlled by
;number of steps).

	Bimf = [initinput(0),initinput(1),initinput(2)]
	Rbs = initinput(3)
	Rmp = initinput(4)
	Vsw = initinput(5)
	dens = initinput(6)
	Xgsminit = initinput(7)
	Ygsminit = initinput(8)
	Zgsminit = initinput(9)
	threshold = initinput(10)

	Xgsm = ptopass(0)
	Ygsm = ptopass(1)
	Zgsm = ptopass(2)

	P0 = ptopass

;extras

	Bms    = [extras(0),extras(1),extras(2)]
	Bearth = [extras(3),extras(4),extras(5)]
	num = extras(6)

;function

	jhatsteps = 0.16
	numsteps = num
	allowbq = threshold
   	coordML = fltarr(3,2.0*numsteps)

;INCREMENT P in the +jhat direction

	NewP = P0

	for i = 0,numsteps-1 do begin

		OldP = NewP
		coordML(*,i) = OldP

		ptopass = OldP

; Following section altered by Robert Fear. The _local_ magnetosheath and magnetospheric magnetic fields at each new point on the X-line
; is calculated using the existing functions, and used when the function 'testing' is called.

		SF = sheathfield(initinput,ptopass)
		Bms_here = SF[0:2] ;FC(sheathfield)
		EF = earthfield(initinput,ptopass, [Rmp/2.0,0.0,Rmp,1.0])  ;FC(earthfield)
		Bearth_here = SF[3]*EF

		extras = [Bms_here,Bearth_here]
		detailsOld = testing(initinput,ptopass,extras) ;FC(testing)
		;tests = [qhat, gradS, jhat, deltaBq, initrecon]

		jhatOld = [detailsOld(6),detailsOld(7), detailsOld(8)]

	;Increment OldP along jhat

		P1 = OldP + jhatOld*jhatsteps

;Bring P1 back to the p'boid

		P2 = lamda(P1, Rmp) ;FC(lamda)

;Calculate details for P2

		ptopass = P2
; Also altered by RCF (see above):
		SF = sheathfield(initinput,ptopass)
		Bms_here = SF[0:2] ;FC(sheathfield)
		EF = earthfield(initinput,ptopass, [Rmp/2.0,0.0,Rmp,1.0])  ;FC(earthfield)
		Bearth_here = SF[3]*EF

		extras = [Bms_here,Bearth_here]
		shear_angle_here = ACOS(dotv1v2(normalise(Bms_here),normalise(Bearth_here)))/!DTOR


		detailsP2 = testing(initinput,ptopass,extras) ;FC(testing)

		jhatP2 = [detailsP2(6),detailsP2(7), detailsP2(8)]

   		deltaBqP2 = abs(detailsP2(9))

	IF anglethreshold THEN BEGIN
		if (shear_angle_here LT allowbq) then NewP = OldP else NewP = P2
	ENDIF ELSE BEGIN
		if (deltabqP2 LT allowbq) then NewP = OldP else NewP = P2
	ENDELSE
	endfor

;j in the -jhat direction

NewP = P0

	for i = numsteps,(2.0*numsteps)-1 do begin

		OldP = NewP
		coordML(*,i) = OldP
; Also altered by RCF (see above):
		ptopass = OldP
		SF = sheathfield(initinput,ptopass)
		Bms_here = SF[0:2] ;FC(sheathfield)
		EF = earthfield(initinput,ptopass, [Rmp/2.0,0.0,Rmp,1.0])  ;FC(earthfield)
		Bearth_here = SF[3]*EF
		extras = [Bms_here,Bearth_here]
		detailsOld = testing(initinput,ptopass,extras) ;FC(testing)
		jhatOld = [detailsOld(6),detailsOld(7), detailsOld(8)]

;Increment OldP along jhat

		P1 = OldP - jhatOld*jhatsteps

;Bring P1 back to the p'boid

		P2 = lamda(P1, Rmp) ;FC(lamda)

;Calculate details for P2

		ptopass = P2
; Also altered by RCF (see above):
		SF = sheathfield(initinput,ptopass)
		Bms_here = SF[0:2] ;FC(sheathfield)
		EF = earthfield(initinput,ptopass, [Rmp/2.0,0.0,Rmp,1.0])  ;FC(earthfield)
		Bearth_here = SF[3]*EF

		extras = [Bms_here,Bearth_here]
		shear_angle_here = ACOS(dotv1v2(normalise(Bms_here),normalise(Bearth_here)))/!DTOR

		detailsP2 = testing(initinput,ptopass,extras) ;FC(testing)

		jhatP2 = [detailsP2(6),detailsP2(7), detailsP2(8)]

   		deltaBqP2 = abs(detailsP2(9))

	IF anglethreshold THEN BEGIN
		if (shear_angle_here LT allowbq) then NewP = OldP else NewP = P2
	ENDIF ELSE BEGIN
		if (deltabqP2 LT allowbq) then NewP = OldP else NewP = P2
	ENDELSE

	;if (deltabqP2 LT allowbq) then NewP = OldP

	;if (deltabqP2 GE allowbq) then NewP = P2

   	endfor

	myline = coordML

;it returns an array containing the coordinates of the merging line to plot on the chart

	return, myline
end