;************************************************************************
The beginnings of an implementation of the Low-Hundhausen Cartesian flux rope model.

At the moment, the prominence current sheet (Aprom) is not included.

When Adip is set to zero (default), observable properties such as

        cavity radius, height, external vertical field, external density at the photosphere 
	and the ratio of density at the flux rope axis to surrounding hydrostatic density (same height)

fully constrain model parameters.


Allowing Adip and eventually Aprom to vary may give some freedom in flux rope twist profile, for example, and possibly also cavity depletion profile.

NOTE -- the temperature profile is increasing with height in this model, which is maybe ok for low 
coronal (below temperature maximum).  However, we have a flag "isothermal" which allows fixing a temperature 
above the bubble (fixed to value at z=ro), which is useful for the LOS integrated quantities.  
If isothermal is set to a value, then the temperature is fixed everywhere (including below the bubble top) to that value.

NOTE -- the model is 2D Cartesian -- it is fit into the spherical coordinates of the FORWARD codes by 
assuming invariance in Phi.  This introduces an error in the magnetic force balance, and there is also 
an error in the graviational description which is gz.  Both should be minimal for small-radius flux rope positioned close to photosphere.
;************************************************************************
;
;
;;Keyword Inputs:
;
;
;               R_OINPUT, X_OINPUT -- these are the radius of the flux rope, 
;			and the percent distance from the axis 
;			to the photosphere (at central meridian, Y=0)
;			both are input in units of solar radii
;			Note that X_O will be negative, with minimum
;			value required to be greater than -1.
;                       DEFAULT R_O = 0.025, X_O = -0.6
;			OUTPUT changed to CGS 
;			Determine ZO,ZETA_O (see below)
;
;		BXEX -- external bipolar field at base external to rope
;			DEFAULT 10 Gauss
;
;		SIGNLAMBDA -- direction of twist
; 			this doesnt affect much except allows
; 			bz/by to be either sign
;			DEFAULT positive (1), OUTPUT UNCHANGED
;	
;		RHOEX -- number density
; 			at base outside flux rope
;			DEFAULT 5d8 (number density)
;			OUTPUT changed to cgs
;
;		RHORAT -- density depletion of cavity
;			DEFAULT 0.5
;
;		ADIP -- external dipole field
;			FREE PARAMETER
;			DEFAULT 0 -- makes problem fully observationally constrained
;			OUTPUT UNCHANGED
;		ISOTHERMAL -- forces temperature to be isothermal for (1) heights
;			above r=ro, with a value = the value at that height.
;			Or if set to a value other than 1, sets the temperature equal to that value everywhere.
;			note -- this screws up the polytropic condition on the model,
;			which creates a increasing temperature with height;  
;			while this is fine for the low corona where the rope sits, 
;			LOS integration would have temperature increasing too much
;
;		HYDRO - forces hydrostatic density falloff -- similar sort of kluge.  Hydrostatic falloff taken to match Lites low.  Value of hydro used to scale external density (hydro=1 is uniform density falloff)
;
;		OUT - this will scale the field strength outside - again a kluge to better compare to Lites Low
			
;               AA,BB,CC,DD,EE,FF - background density, pressure parameters,
;                       for hydrostatic equilibrium OUTPUT as above but
;                       aa, cc, ee  multiplied by 1d8
;                       see below for DEFAULT
;
;
;               VELIMPOSE - impose a velocity of magnitude VELIMPOSE 
;               directed along the field
;               will overwrite any velocity field already loaded into the cube
;               if nonzero
;                       DEFAULT 0.d0