function NuclearChargeDensity()
{
	this.Const = new Constants();
	this.cf = new CommonFun();
	this.gs = new GaussQuadrature();
	
	this.Debugging = false;
	this.Show = true;
	this.TypeAbbrevation = '';
    this.ShortReference = '';
    this.Reference = '';
	this.Z = -1;
	this.A = -1;
	this.Type = -1;
	this.IndexRef = -1;

	this.R2 = 0;
	this.RP = 0;
	this.c = 0;
	this.a = 0;
	this.z = 0;
	this.alpha = 0;
	this.w = 0;
	this.Qmin = 0;
	this.Qmax = 1;
	this.Ri = [];
	this.Qi = [];
	this.Ai = [];
	this.rho0 = 0;

this.defineRho0 = function()
{
	this.gs.setGaussQuadrature(20);	
	this.rho0 = 1.0;
    var 
		Integral = 0,
		ig,
        r,rho,
		rmax = 3.0*this.R2;
    for(ig=0;ig < this.gs.NGauss; ig++)
    {
        r = rmax*0.5*(1*this.gs.xGauss[ig]+1.0);
		rho = this.getDensity(r);
        Integral = 1*Integral + 1*r*r*this.gs.wGauss[ig]*rho;
    }
    this.rho0 = this.Z/(4*Math.PI*Integral*0.5*rmax);
}
this.getAnaliticExpression = function()
{
	var buffer = '',i;
    switch(this.Type)
    {
        default: break;
        case 0: // Harmonic oscilator
			buffer += "Harmonic oscilator model:\n"+
				"\\rho(r) = \\rho_0(1+\\alpha r^2/a^2) exp(-r^2/a^2),\n"+
	        	"where a = "+this.cf.stringInFormat(this.a,4)+" fm, "+
	        	"\\alpha = "+this.cf.stringInFormat(this.alpha,4)+", "+
	        	"\\rho_0 = "+this.cf.stringInFormat(this.rho0,4)+" fm^{-3}";
            break;
        case 1: // Modified Harmonic oscilator
			buffer += "Modified harmonic oscilator model:\n"+
				"\\rho(r) = \\rho_0(1+\\alpha r^2/a^2) exp(-r^2/a^2),\n"+
	        	"where a = "+this.cf.stringInFormat(this.a,4)+" fm, "+
	        	"\\alpha = "+this.cf.stringInFormat(this.alpha,4)+
	        	"\\rho0 = "+this.cf.stringInFormat(this.rho0,4)+" fm^{-3}";
            break;
        case 5: // 2 parameters Fermi model
			buffer += "2 parameters Fermi model:\n"+
				"\\rho(r) = \\rho_0/(1+exp((r-c)/z)),\n"+
	        	"where c = "+this.cf.stringInFormat(this.c,4)+" fm, "+
	        	"z = "+this.cf.stringInFormat(this.z,4)+" fm, "+
	        	"\\rho0 = "+this.cf.stringInFormat(this.rho0,4)+" fm^{-3}";
            break;
        case 6: // 3 parameters Fermi model
			buffer += "3 parameters Fermi model:\n"+
				"\\rho(r) = \\rho_0(1+w r^2/c^2)/(1+exp((r-c)/z)),\n"+
	        	"where c = "+this.cf.stringInFormat(this.c,4)+" fm, "+
	        	"z = "+this.cf.stringInFormat(this.z,4)+" fm, "+
	        	"w = "+this.cf.stringInFormat(this.w,4)+", "+
	        	"\\rho0 = "+this.cf.stringInFormat(this.rho0,4)+" fm^{-3}";
            break;
        case 7: // 3 parameters Gauss model
			buffer += "3 parameters Gaussian model:\n"+
				"\\rho(r) = \\rho_0(1+w r^2/c^2)/(1+exp((r^2-c^2)/z^2)),\n"+
	        	"where c = "+this.cf.stringInFormat(this.c,4)+" fm, "+
	        	"z = "+this.cf.stringInFormat(this.z,4)+" fm, "+
	        	"w = "+this.cf.stringInFormat(this.w,4)+", "+
	        	"\\rho0 = "+this.cf.stringInFormat(this.rho0,4)+" fm^{-3}";
            break;
        case 3: // Fourier-Bessel
			buffer += "Fourier-Bessel expansion:\n"+
				"\\rho(r) = \\sum_\\nu a_\\nu j_0(\\nu\\pi r/R) for r < R,\n"+
				"\\rho(r) = 0.0 for r > R,\n"+
	        	"here R = "+this.cf.stringInFormat(this.RP,4)+" fm, and \n";
	        for(i=0;i<15;i+=3)
				buffer += "a_{"+(i+1)+"}="+this.Ai[i]+", a_{"+(i+2)+"}="+this.Ai[i+1]+", a_{"+(i+3)+"}="+this.Ai[i+2]+",\n";
			buffer += "a_{16}="+this.Ai[15]+", a_{17}="+this.Ai[16];
			break;
        case 4: // Gauss sum
			buffer += "Sum-of-Gaussians parameters:\n"+
				"\\rho(r) = \\sum_i A_i [exp(-(r-R_i)^2/\\gamma^2) + exp(-(r+R_i)^2/\\gamma^2)],\n"+
	        	"where A_i = Z Q_i/[2\\pi^{3/2}\\gamma^3 (1+2R_i^2/\\gamma^2)], and \n"+
	        	"\\gamma = \\sqrt{2/3} R_G, and the rms radius of the Gaussians R_G="+this.cf.stringInFormat(this.RP,4)+" fm,\n";
	        for(i=1;i<=12;i++) buffer += "R_{"+i+"}="+this.Ri[i-1]+" Q_{"+i+"}="+this.Qi[i-1]+"\n";
			break;
        case 8: // Uniform Gauss model
            break;
    }
    buffer += "\n\n";
    return buffer;
}
this.getDensity = function(r)
{
    var
		i,ra,r1,r2,c2,z2,
        Ci,gamma2,
        rho = 0;
    switch(this.Type)
    {
        default: break;
        case 0: // Harmonic oscilator
        case 1: // Modified Harmonic oscilator
            ra = r*r/(this.a*this.a);
            if(ra>1.e100) rho = 0;
            else rho = this.rho0*(1+1*this.alpha*ra)*Math.exp(-ra);
            break;
        case 5: // 2 parameters Fermi model
            rho = this.rho0/(1+1*Math.exp((r-this.c)/this.z));
            break;
        case 6: // 3 parameters Fermi model
            rho = this.rho0*(1+1*this.w*r*r/(this.c*this.c))/(1+1*Math.exp((r-this.c)/this.z));
            break;
        case 7: // 3 parameters Gauss model
            r2 = r*r; c2 = this.c*this.c; z2 = this.z*this.z;
            rho = this.rho0*(1+1*this.w*r2/c2)/(1+1*Math.exp((r2-c2)/z2));
            break;
        case 2: // Model independent evaluation
        case 3: // Fourier-Bessel
            rho = 0.0;
            if(r<this.RP)
            {
                for(i=1;i<=17;i++)
                {
                    r1 = i*Math.PI*r/this.RP;
                    if(r1<1.e-6) r2 = 1.0; else r2 = Math.sin(r1)/r1;
                    rho = 1*rho + 1*this.Ai[i-1]*r2;
                }
                rho *= this.rho0;
            }
            break;
        case 4: // Gauss sum
            rho = 0.0;
            gamma2 = this.RP*this.RP*0.666666666667;
            c2 = this.Z*0.5*this.Const.M_1_SQRT_PI*this.Const.M_1_PI/(gamma2*Math.sqrt(gamma2));
            for(i=0;i<12;i++)
            {
                Ci = this.Qi[i]*c2/(1.0+2.0*this.Ri[i]*this.Ri[i]/gamma2);
                r1 = (r-this.Ri[i])*(r-this.Ri[i])/gamma2;
                r2 = (1*r+1*this.Ri[i])*(1*r+1*this.Ri[i])/gamma2;
                rho = 1*rho + 1*Ci*(Math.exp(-r1) + Math.exp(-r2));
            }
            rho *= this.rho0;
            break;
        case 8: // Uniform Gauss model
            break;
    }
    return rho;
}
this.setNucleus = function(zz,aa)
{
    this.Z = zz;
    this.A = aa;
}
this.setFourierBesselParameters = function(r2,rp,ai)
{
    var i;
	this.R2 = r2; 
	this.RP = rp;
    this.Qmin = 0.0;
    this.Qmax = 0.0;
    this.Ai = [];
    switch(this.Type)
    {
        default: return;
        case 3: // Fourier-Bessel
            for(i=0;i<17;i++) this.Ai[i] = ai[i];
            break;
    }
    this.defineRho0();
}
this.setGaussParameters = function(r2,rp,ri,qi)
{
    var i;
	this.R2 = r2; this.RP = rp;
    this.Qmin = 0.0;
    this.Qmax = 0.0;
    this.Ri = [];
    this.Qi = [];
    switch(this.Type)
    {
        default: return;
        case 4: // Gauss sum
            for(i=0;i<12;i++)
            {
                this.Ri[i] = ri[i];
                this.Qi[i] = qi[i];
            }
            break;
    }
    this.defineRho0();
}
this.setUniformParameters = function(r2,c_a,z_alpha,ww,qmin,qmax)
{
    this.R2 = r2;
    this.Qmin = qmin;
    this.Qmax = qmax;
    switch(this.Type)
    {
        default: return;
        case 0: // Harmonic oscilator
        case 1: // Modified Harmonic oscilator
            this.a = c_a;
            this.alpha = z_alpha;
            break;
        case 5: // 2 parameters Fermi model
            this.c = c_a;
            this.z = z_alpha;
            break;
        case 6: // 3 parameters Fermi model
        case 7: // 3 parameters Gauss model
            this.c = c_a;
            this.z = z_alpha;
            this.w = ww;
            break;
        case 2: // Model independent evaluation
        case 8: // Uniform Gauss model
            break;
    }
    this.defineRho0();
}
this.setReference = function(ref)
{
	this.ShortReference = ref;
}
this.setShortReference = function(ref)
{
	this.ShortReference = ref;
}
this.setType = function(type)
{
	type = type.trim();
	this.TypeAbbrevation = type;
    
    if(type == "HO")       this.Type = 0;
    else if(type == "MHO") this.Type = 1;
    else if(type == "MI")  this.Type = 2;
    else if(type == "FB")  this.Type = 3;
    else if(type == "GS")  this.Type = 4;
    else if(type == "2pF") this.Type = 5;
    else if(type == "3pF") this.Type = 6;
    else if(type == "3pG") this.Type = 7;
    else if(type == "UG")  this.Type = 8;
    else this.Type = -1;
}
}