/*	Objective - expression for :  /%
1)			ETods(POIo,t)  =  E0ods(POIo,t)  +  EIods(POIo,t)  

/*	Gaussian coordinates - see also 04_33 Equation (7), see [explanation, worry] below  /%
2)			E0ods(POIo,t)  =  Q(PART)*Rpcs(POIp)^(-2)

/*	EIods(POIo,t,2nd stage) - starting with 04_34 differentiable result   /%
3)	EIods(POIo,t,2nd stage) 
	=	+	K_1st
		+	β*Rocs(POIo)^2
			*{	+ 21/2						*β^3*Q(PART)*Rocs(POIo)^2	*Rpcs(POIo(t),t)^(-6)*∫[∂(Aθpc),0 to Aθpc(POIo(t),t): sin(Aθpc(POIo(t),t))^3	*cos(Aθpc(POIo(t),t)) ]
					- 2	*λ(Vons(PART))	*β	 *Q(PART)/Rocs(POIo)		*Rpcs(POIo(t),t)^(-3)*∫[∂(Aθpc),0 to Aθpc(POIo(t),t): sin(Aθpc(POIo(t),t))		*cos(Aθpc(POIo(t),t)) ]
			 }
		+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}

/*_file_insert_path "$d_Lucas""relativistic factor, intermediate symbols.txt"

/*	using   /%
2687:(mathH)  ∫[∂(Aθpc),0.to.Aθpc(POIo(t),t=0): cos(Aθpc(POIo(t),t))*sin(Aθpc(POIo(t),t)))^3]  = sin(Aθpc(POIo(t),t=0))^4/4
2685:(mathH)  ∫[∂(Aθpc),0.to.Aθpc(POIo(t),t=0): cos(Aθpc(POIo(t),t))*sin(Aθpc(POIo(t),t)))^1]  = sin(Aθpc(POIo(t),t=0))^2/2 

4)		EIods(POIo,t=0,2nd stage)  
		=	+	K_1st 
			+	β*Rocs(POIo)^2
				*{	+ 21/2						*β^3*Q(PART)*Rocs(POIo)^2	*Rpcs(POIo(t),t)^(-6)*sin(Aθpc(POIo(t),t=0))^4/4
						- 2	*λ(Vons(PART))	*β	 *Q(PART)/Rocs(POIo)		*Rpcs(POIo(t),t)^(-3)*sin(Aθpc(POIo(t),t=0))^2/2 
				 }
			+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}
5)		EIods(POIo,t=0,2nd stage)  
		=	+	K_1st 
			+	β*Rocs(POIo)^2	*21/2*β^3	*Q(PART)*Rocs(POIo)^2	*Rpcs(POIo(t),t)^(-6)*sin(Aθpc(POIo(t),t=0))^4/4
			-	β*Rocs(POIo)^2	*2	  *β		*Q(PART)/Rocs(POIo)		*Rpcs(POIo(t),t)^(-3)*sin(Aθpc(POIo(t),t=0))^2/2 *λ(Vons(PART))	
			+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}
/*	collect factors, rearrange using  /%
1108:(mathH)/*	differentiable form  /%
	K_1st 
	=	+	Q(PART)
			*(							3/2	*β^2*Rocs(POIo)^3	*Rpcs(POIo(t),t)^(-5)*sin(Aθpc(POIo(t),t))^2 
				-	λ(Vons(PART))			 						*Rpcs(POIo(t),t)^(-2)
			 )

6)		EIods(POIo,t=0,2nd stage) 
		=	+	Q(PART)
				*(							3/2	*β^2*Rocs(POIo)^3	*sin(Aθpc(POIo(t),t=0))^2	*Rpcs(POIo(t),t=0)^(-5)
					+ 						21/8	*β^4*Rocs(POIo)^4	*sin(Aθpc(POIo(t),t=0))^4	*Rpcs(POIo(t),t=0)^(-6)
					-	λ(Vons(PART))	*1	  														 		*Rpcs(POIo(t),t=0)^(-2)
					-	λ(Vons(PART))	*1		*β^2*Rocs(POIo)	*sin(Aθpc(POIo(t),t=0))^2	*Rpcs(POIo(t),t=0)^(-3)
				 )
			+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}

/*	what Lucas seems to have done here is to use :  /%
1205:(mathH)  E0ods(POIo,t) =  E0pds(POIo(t),t) = Q(PART)/Rpcs(POIo(t),t)^2
(mathL)/*	HIGHLY restrictive conditions : t=0, RFo(t=0) = Rfp(t=0), Rocs(POIo) = Rpcs(POIo(t),t=0)  /%
	Rocs(POIo) = Rpcs(POIo(t),t=0)
(endMath)
/*	|--???-->  
	02Sep2019  Note that I find that the assumption {Rocs(POIo) = Rpcs(POIo(t),t=0)} is a HUGE restriction, 
	and means to me that the result is not at all general (but neither are standard theories).
	<--???--|
/*	re-express (6) using Rocs(POIo) = Rpcs(POIo(t),t=0)  /%
7)		EIods(POIo,t=0,2nd stage) 
		=	+						3/2	*β^2*Q(PART)*Rpcs(POIo(t),t=0)^(-2)	*sin(Aθpc(POIo(t),t=0))^2
			+						21/8	*β^4*Q(PART)*Rpcs(POIo(t),t=0)^(-2)	*sin(Aθpc(POIo(t),t=0))^4
			-	λ(Vons(PART))	*1			 *Q(PART)*Rpcs(POIo(t),t=0)^(-2)
			-	λ(Vons(PART))	*1		*β^2*Q(PART)*Rpcs(POIo(t),t=0)^(-2)	*sin(Aθpc(POIo(t),t=0))^2
			+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}
/*	rearranging  /%
8)		EIods(POIo,t=0,2nd stage) 
		=	+	Q(PART)*Rpcs(POIo(t),t=0)^(-2)
				*(	+						3/2	*β^2*sin(Aθpc(POIo(t),t=0))^2
					+						21/8	*β^4*sin(Aθpc(POIo(t),t=0))^4
					-	λ(Vons(PART))	*1 
					-	λ(Vons(PART))	*1		*β^2*sin(Aθpc(POIo(t),t=0))^2
				 )
			+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}

/*	now using  /%
1174:(mathH)  E0pds(POIp) = Q(PART)/Rpcs(POIp)^2
/*	substitute (2) into (8), and (8) into (1), drop f_sphereCapSurf expression  /%
9)		EIods(POIo,t=0,2nd stage)  
		=	+ E0ods(POIo,t)						*3/2	*β^2*sin(Aθpc(POIo(t),t=0))^2
			- E0ods(POIo,t)	*λ(Vons(PART))	*1
			+ E0ods(POIo,t)						*21/8	*β^4*sin(Aθpc(POIo(t),t=0))^4
			- E0ods(POIo,t)	*λ(Vons(PART))	*1		*β^2*sin(Aθpc(POIo(t),t=0))^2	
/*	re-express (9)  /%
10)	EIods(POIo,t=0,2nd stage)  
		=  + E0ods(POIo,t)					*{		 3/2*β^2*sin(Aθpc(POIo(t),t=0))^2 + 21/8*β^4*sin(Aθpc(POIo(t),t=0))^4 }
			- E0ods(POIo,t)*λ(Vons(PART))	*{ 1 + 	  β^2*sin(Aθpc(POIo(t),t=0))^2 }
/*	For total E, ETods(POIo,t)  /%
11)	ETods(POIo,t=0,2nd stage)  
		=  E0ods(POIo,t)  +  EIods(POIo,t) 
		=  + E0ods(POIo,t)					*{ 1 + 3/2*β^2*sin(Aθpc(POIo(t),t=0))^2 + 21/8*β^4*sin(Aθpc(POIo(t),t=0))^4 }
			- E0ods(POIo,t)*λ(Vons(PART))	*{ 1 + 	  β^2*sin(Aθpc(POIo(t),t=0))^2 }


/*++++++++++++++++++++++++++++++++++++++
/*add_eqn  03Sep2019 NOT the same thing!!!
   04_35
   F therefore E balance - more iterations 
	all measures at t=0 
 
/$L	Ei(r)|t=0  
		= + E0(rs)     *(    3/2*β^2*sin(θ)^2 + 15/8*β^4*sin(θ)^4 )  
		  - E0(rs)*λ(v)*(1 + 3/2*β^2*sin(θ)^2)  

/$H	Ei(rs,vs)  all measures at t=0  
		= + E0(rs)     *(    3/2*β^2*sin(θ´)^2 + 15/8*β^4*sin(θ´)^4 )  
		  - E0(rs)*λ(v)*(1 + 3/2*β^2*sin(θ´)^2)  

/*	OK other than Lucas EI/ET mixup, Oo/Op limit used by Lucas p72h0.0? probably OK, check later
    E0 is added to get ET in 4-36.
    As part of the iterative approach, Lucas dropped the last expression with Ei(r).  
    This is in need of [explanation, clarification] - i.e. this is the first-step approximation only. 
/* Both [induced, total] forms are shown below  
/%H

/*	remove "=0" from t to get : 
(mathL)/*	differentiable form  EIods(POIo,t,2nd stage)  /%
	EIods(POIo,t,2nd stage) 
	=	+	Q(PART)
				*(	+						3/2	*β^2*Rocs(POIo)^3	*sin(Aθpc(POIo(t),t))^2	*Rpcs(POIo(t),t)^(-5)
					+ 						21/8	*β^4*Rocs(POIo)^4	*sin(Aθpc(POIo(t),t))^4	*Rpcs(POIo(t),t)^(-6)
					-	λ(Vons(PART))	*1	  															*Rpcs(POIo(t),t)^(-2)
					-	λ(Vons(PART))	*1		*β^2*Rocs(POIo)	*sin(Aθpc(POIo(t),t))^2	*Rpcs(POIo(t),t)^(-3)
				 )
		+ f_sphereCapSurf{f_sphereCapSurf(EIods(POIo,t,0th stage))]}
(endMath)

(mathL)/*	HIGHLY restrictive conditions : t=0, RFo(t=0) = Rfp(t=0), Rocs(POIo) = Rpcs(POIo(t),t=0)  /%
	EIods(POIo,t=0,2nd stage) 
	=		E0ods(POIo,t)					*{		3/2*β^2*sin(Aθpc(POIo(t),t=0))^2 + 21/8*β^4*sin(Aθpc(POIo(t),t=0))^4}
		-	E0ods(POIo,t)*λ(Vons(PART))*{1 +		 β^2*sin(Aθpc(POIo(t),t=0))^2}
(endMath)

(mathL)/*	HIGHLY restrictive conditions : t=0, RFo(t=0) = Rfp(t=0), Rocs(POIo) = Rpcs(POIo(t),t=0)  /%
	ETods(POIo,t=0,2nd stage)
	=		E0ods(POIo,t)					*{1 + 3/2*β^2*sin(Aθpc(POIo(t),t=0))^2 + 21/8*β^4*sin(Aθpc(POIo(t),t=0))^4} 
		-	E0ods(POIo,t)*λ(Vons(PART))*{1 +		 β^2*sin(Aθpc(POIo(t),t=0))^2}
(endMath)

/*_file_insert_path "$d_Lucas""relativistic factor, restrictive conditions.txt"

Note that Lucas (4-35) is actually for EIods(POIo,t=0), not ETods(POIo,t).
/*	03Sep2019 the numerical factors differ from Lucas - hopefully that will self-fix with further iterations?