def theta_generate(system, plot = False, realonly = False, useapproximations = False):
        global theta_2,theta_3, cos_theta_2, cos_theta_3, tan_theta_2
        (n_1,n_2,n_3) = system.indices
        
        theta_1 = system.theta_1
        freq = system.freq


        cos_theta_2 = (n_2)**-1  * array(n_2**2-((n_1)*sin(theta_1))**2,dtype = complex)
        cos_theta_3 = (n_3)**-1  * array(n_3**2-((n_1)*sin(theta_1))**2,dtype = complex)
        tan_theta_2 = sqrt((1-cos_theta_2**2)/cos_theta_2)

            
        theta_2 = numpy.arcsin(sin(theta_1) * n_1/ n_2,dtype = complex)
        theta_3 = numpy.arcsin(sin(theta_1)* n_1 / n_3,dtype = complex)

        
        
        if plot == True:

                
                py.plot(freq,numpy.real(cos_theta_2[IR]),'r')
                py.plot(freq,numpy.imag(cos_theta_2[IR]),'b')
                py.plot(freq,numpy.real(cos_theta_3[IR]),'g')
                py.plot(freq,numpy.imag(cos_theta_3[IR]),'k')
                py.title("cos Angles of incidence")
                py.legend(('real2','imag2','real3','imag3'))
                py.show()
        
        system.angles = (theta_2,theta_3, cos_theta_2, cos_theta_3, tan_theta_2)
        
        return 0

 def __arc_segment(self,r,a,b):
   """
     Returns length of circle segment in a half-infinite slit.
       r   ... radius of circle which is centered at the origin
       a,b ... borders of slit containing (x,y) with x>a, b>y>0 
   """  
   #print r, a, b
   # scalar version
   if np.isscalar(r):
     if r<a:              # circle outside slit
       return 0;   
     elif r**2>a**2+b**2: # circle does not intersect left border at x=a
       return r*np.arcsin(b/r);
     else:                # circle only intersects left border at x=a
       return r*np.arccos(a/r);
 
   # parallel version
   else:
     r     = np.atleast_1d(r); 
     arc   = np.zeros(r.size);
     index = r>=a;              # select radii which intersect half-slit
     y2_a  = r[index]**2-a**2;  # intersection of circle with left border at x=a
     ymax  = np.where( y2_a > b**2, b, np.sqrt(y2_a) ); 
                                # largest y value on circle segment in half-slit
     arc[index] = r[index] * np.arcsin( ymax / r[index] );
                                # calculate arc length
     # DEBUG
     if self.verbosity>2:
       arc2 = [ self.__arc_segment(ir,a,b) for ir in r ];
       assert(np.allclose( arc, arc2 ));
     return arc;
 
   raise RuntimeError("Should never reach this point");

 def getCurvePoints(self,P0,P1,P2,R,n=20):
     """Generates a curve. P0,P1 and P2 define a corner and R defines the radius of the
     tangential circle element. n gives the number of point to approximate the curve.
     P1 is the corner itself and P0 and P2 give the tangents.
     """
     P0=np.array(P0)
     P1=np.array(P1)
     P2=np.array(P2)
     o1=(P0-P1)/np.sqrt(np.dot(P0-P1,P0-P1))
     o2=(P2-P1)/np.sqrt(np.dot(P2-P1,P2-P1))
     if np.arcsin(np.cross(o1,o2)) > 0:
         a=1.
         b=-1.
     else:
         a=-1.
         b=1.
     
     v1=R*np.dot(np.array([[0.,b],[a,0.]]),o1)
     v2=R*np.dot(np.array([[0.,a],[b,0.]]),o2)
     dv=v2-v1
     a=np.array([[o1[0],-o2[0]],[o1[1],-o2[1]]])
     b=dv
     x=np.linalg.solve(a,b)
     circleCenter= P1+x[0]*o1+v1
     angle = np.arcsin(np.cross(v2/R,v1/R))
     points=[]
     for i in range(n+1):
         x=-i*angle/n
         rot = np.array([[np.cos(x),-np.sin(x)],
                          [np.sin(x),np.cos(x)]])
         points.append(circleCenter+np.dot(rot,-v1))
     return points

    def __remove_particles_randomly(self):

        grid_pairs = np.array([self.grid[0].ravel(), self.grid[1].ravel()]).T

        particles = grid_pairs[
            grid_pairs[:, 1]
            > (
                self.sin_amp * np.sin(self.sin_freq * grid_pairs[:, 0] + np.arcsin(1.0))
                + self.sin_amp
                + self.particle_diameter
            )
        ]
        self.driver = grid_pairs[
            grid_pairs[:, 1]
            <= (
                self.sin_amp * np.sin(self.sin_freq * grid_pairs[:, 0] + np.arcsin(1.0))
                + self.sin_amp
                + self.particle_diameter
            )
        ]

        self.x = particles[:, 0]

        number_of_particles_to_remove = self.__compute_number_of_particles_to_remove()

        np.random.shuffle(particles)

        self.x = particles[:-number_of_particles_to_remove, 0]
        self.y = particles[:-number_of_particles_to_remove, 1]

def convert_galactic_equa(l, b):
   '''l and b in degrees'''
   r = pi / 180.0
   dec = arcsin(cos(b * r) * sin((l - 33) * r) * sin(62.6 * r) + sin(b * r) * cos(62.6 * r))
   ra  = 282.25 + arcsin((cos(b * r) * sin((l - 33) * r) * cos(62.6 * r) - sin(b * r) * sin(62.6 * r)) / cos(dec)) / r
   dec = dec / r
   return ra, dec

  def test_cartesianToSpherical(self):
    """
    Check correct working of the Cartesian to Spherical coordinates transformation.
    """
    r, phi, theta = cartesianToSpherical(array([1,0,0]), array([0,1,0]), array([0,0,1]))
    assert_array_almost_equal(r,array([1.0,1.0,1.0]))
    assert_array_almost_equal(phi,array([0.0,pi/2.0,0.0]))
    assert_array_almost_equal(theta,array([0.0,0.0,pi/2.0]))

    r, phi, theta = cartesianToSpherical(1.0, 1.0, 0.0)
    assert_almost_equal(r, sqrt(2.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, 0.0)

    r, phi, theta = cartesianToSpherical(1.0, 1.0, 1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, arcsin(1.0/sqrt(3.0)))

    r, phi, theta = cartesianToSpherical(1.0, 1.0, -1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0)
    assert_almost_equal(theta, -arcsin(1.0/sqrt(3.0)))

    r, phi, theta = cartesianToSpherical(-1.0, -1.0, -1.0)
    assert_almost_equal(r, sqrt(3.0))
    assert_almost_equal(phi, pi/4.0+pi)
    assert_almost_equal(theta, -arcsin(1.0/sqrt(3.0)))

    assert_raises(Exception, cartesianToSpherical, 0.0, 0.0, 0.0)
    assert_raises(Exception, cartesianToSpherical, array([1,0,0,0]), array([0,1,0,0]), array([0,0,0,1]))

def solve_acos_plus_bsin_plus_c(a, b, c):
    """
    Solves a*cos(x) + b*sin(x) + c = 0 for x in [-pi/2, pi/2].

    The method used is to change variable by setting t = sin(x).

    Args:
        a, b, c: coefficients of the equation

    Returns:
        A solution x if it exists, None in the other cases.
    """
    out = None
    if a * a + b * b - c * c > 0:
        poly = np.poly1d([a * a + b * b, 2 * b * c, c * c - a * a])
        if not any([-1 <= t <= 1 for t in poly.r]):
            print("no solutions between -1 and 1: ", poly.r)
        elif all([-1 <= t <= 1 for t in poly.r]):
            i = np.argmin(np.abs(np.array([a * np.sqrt(1 - t * t) + b * t + c for t in poly.r])))
            out = np.arcsin(poly.r[i])
        else:
            if -1 <= poly.r[0] <= 1:
                out = np.arcsin(poly.r[0])
            else:
                np.arcsin(poly.r[1])

    # check that the solution is in [-pi/2, pi/2] and satisfies the equation
    assert -np.pi / 2 < out < np.pi / 2
    np.testing.assert_allclose(a * np.cos(out) + b * np.sin(out) + c, 0, atol=1e-12)
    return out

def tthToRad(twoTheta, unit, wavelength=None, directDist=None):
    """
    Convert a two theta angle from original `unit` to radian.

    `directDist = ai.getFit2D()["directDist"]`
    """
    if isinstance(twoTheta, numpy.ndarray):
        pass
    elif isinstance(twoTheta, collections.Iterable):
        twoTheta = numpy.array(twoTheta)

    if unit == units.TTH_RAD:
        return twoTheta
    elif unit == units.TTH_DEG:
        return numpy.deg2rad(twoTheta)
    elif unit == units.Q_A:
        if wavelength is None:
            raise AttributeError("wavelength have to be specified")
        return numpy.arcsin((twoTheta * wavelength) / (4.e-10 * numpy.pi)) * 2.0
    elif unit == units.Q_NM:
        if wavelength is None:
            raise AttributeError("wavelength have to be specified")
        return numpy.arcsin((twoTheta * wavelength) / (4.e-9 * numpy.pi)) * 2.0
    elif unit == units.R_MM:
        if directDist is None:
            raise AttributeError("directDist have to be specified")
        # GF: correct formula?
        return numpy.arctan(twoTheta / directDist)
    elif unit == units.R_M:
        if directDist is None:
            raise AttributeError("directDist have to be specified")
        # GF: correct formula?
        return numpy.arctan(twoTheta / (directDist * 0.001))
    else:
        raise ValueError("Converting from 2th to unit %s is not supported", unit)

    def reconstruct_angle(self, event, offsets=None):
        """Reconstruct angles from a single event"""

        c = 3.00e+8

        if offsets is not None:
            self._correct_offsets(event, offsets)

        dt1 = event['t1'] - event['t3']
        dt2 = event['t1'] - event['t4']

        station = self.cluster.stations[event['station_id']]
        r1, phi1 = station.calc_r_and_phi_for_detectors(1, 3)
        r2, phi2 = station.calc_r_and_phi_for_detectors(1, 4)

        phi = arctan2((dt2 * r1 * cos(phi1) - dt1 * r2 * cos(phi2)),
                      (dt2 * r1 * sin(phi1) - dt1 * r2 * sin(phi2)) * -1)
        theta1 = arcsin(c * dt1 * 1e-9 / (r1 * cos(phi - phi1)))
        theta2 = arcsin(c * dt2 * 1e-9 / (r2 * cos(phi - phi2)))

        e1 = sqrt(self.rel_theta1_errorsq(theta1, phi, phi1, phi2, r1, r2))
        e2 = sqrt(self.rel_theta2_errorsq(theta2, phi, phi1, phi2, r1, r2))

        theta_wgt = (1 / e1 * theta1 + 1 / e2 * theta2) / (1 / e1 + 1 / e2)

        if theta_wgt < 0:
            theta_wgt *= -1
            phi += pi
            phi = (phi + pi) % (2 * pi) - pi

        return theta_wgt, phi

def tpi2k(id, p, dI):
    p2brho = 10./2.99792458
    deg = pi/180.0
    if id not in table.keys():
        print( "type identifier is unknown")
        print( "known types:")
        print( list(table.keys())[:10])
        print( list(table.keys())[10:20])
        print( list(table.keys())[20:30])
        print( list(table.keys())[35:])
        return None
    magnet_type = table[id]["type"]
    a = table[id]["A"]
    el = table[id]["EL"]
    km = table[id]["KM"]
    dk0 = p3(a, dI ) * km/(p*p2brho)
    dk = 0.
    if magnet_type == "S" or magnet_type == "Q":
        #case kQ:
        dk = dk0
    elif magnet_type == "B":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi/deg
    elif magnet_type ==  "C":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi*1.e3
    elif magnet_type == "T":
        phi = 2.0*np.arcsin(0.5*el * dk0)
        dk = phi*1.e3
    return dk

def test_rphi_to_dl_2d():
    #This is a tangent point
    r,phi= 6., numpy.arccos(0.75)
    d,l= bovy_coords.rphi_to_dl_2d(r,phi,degree=False,ro=8.,phio=0.)
    l= numpy.arcsin(0.75)
    d= 6./numpy.tan(l)
    assert numpy.fabs(d-6./numpy.tan(numpy.arcsin(0.75))) < 10.**-10., 'dl_to_rphi_2d conversion did not work as expected'
    assert numpy.fabs(l-numpy.arcsin(0.75)) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point
    r,phi= 2., 55.
    d,l= bovy_coords.rphi_to_dl_2d(r,phi,degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=10.)
    assert numpy.fabs(d-2.) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.fabs(l-45.) < 10.**-10., 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point, for arrays
    r,phi= 2., 45.
    os= numpy.ones(2)
    d,l= bovy_coords.rphi_to_dl_2d(os*r,os*phi,
                                   degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=0.)
    assert numpy.all(numpy.fabs(d-2.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.all(numpy.fabs(l-45.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    #This is another point, for lists, which for some reason I support
    r,phi= 2., 45.
    d,l= bovy_coords.rphi_to_dl_2d([r,r],[phi,phi],
                                   degree=True,ro=2.*numpy.sqrt(2.),
                                   phio=0.)
    d= numpy.array(d)
    l= numpy.array(l)
    assert numpy.all(numpy.fabs(d-2.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    assert numpy.all(numpy.fabs(l-45.) < 10.**-10.), 'rphi_to_dl_2d conversion did not work as expected'
    return None

 def _cal_center(p0,p1,p2):
     #p0 is of type up-down 
     origin=(p1+p2)/2
     y_v=f3(np.zeros(3),p1-origin)
     x_v=f3(np.zeros(3),p0-origin)
     z_v=np.cross(x_v,y_v)
     T=f1(x0_v,y0_v,z0_v,x_v,y_v,z_v)
     #print T
     r=f2(p1,p2)/2*np.tan(np.pi/2-self.open_angle/2.)
     phi=0.
     L1=f2(p1,p2)/2./np.tan(self.open_angle/2.)
     L2=f2(p0,origin)
     #look at document#1 in binder for detail
     a,b,c=1+np.tan(self.theta_top_down)**2,2*L1/L2*np.tan(self.theta_top_down),(L1/L2)**2-1
     sin_list=[(b+(b**2-4*a*c)**0.5)/2./a,(b-(b**2-4*a*c)**0.5)/2./a]
     theta=0
     if (sin_list[0]<1)&(sin_list[0]>0):
         if self.switch==False:
             theta=np.pi/2+np.arcsin(sin_list[0])
         elif self.switch==True:
             theta=np.pi/2-np.arcsin(sin_list[0])
     else:
         if self.switch==False:
             theta=np.pi/2+np.arcsin(sin_list[1])
         elif self.switch==True:
             theta=np.pi/2-np.arcsin(sin_list[1])
         
     center_point_new=np.array([r*np.cos(phi)*np.sin(theta),r*np.sin(phi)*np.sin(theta),r*np.cos(theta)])
     center_point_org=np.dot(inv(T),center_point_new)+origin
     return center_point_org

def R2axis(M):

    m00 = M[0, 0]
    m01 = M[0, 1]
    m02 = M[0, 2]
    m10 = M[1, 0]
    m11 = M[1, 1]
    m12 = M[1, 2]
    m20 = M[2, 0]
    m21 = M[2, 1]
    m22 = M[2, 2]
    # symmetric matrix K
    K = np.array(
        [
            [m00 - m11 - m22, 0.0, 0.0, 0.0],
            [m01 + m10, m11 - m00 - m22, 0.0, 0.0],
            [m02 + m20, m12 + m21, m22 - m00 - m11, 0.0],
            [m21 - m12, m02 - m20, m10 - m01, m00 + m11 + m22],
        ]
    )
    K /= 3.0
    # quaternion is eigenvector of K that corresponds to largest eigenvalue
    w, V = np.linalg.eigh(K)
    q = V[[3, 0, 1, 2], np.argmax(w)]

    # quaternion to Axis
    nrm = np.linalg.norm(q[1:4])
    if nrm <= 0.0:
        a = np.array([0.0, 0.0, 0.0])
    elif q[0] < 0.0:
        a = (np.pi - np.arcsin(nrm)) * 2 * q[1:4] / nrm
    else:
        a = np.arcsin(nrm) * 2 * q[1:4] / nrm

    return a

def funDR1R2(x, v1, v2, v3, z1, z2):
# def funDR1R2(x, v1, v2, v3, z1, z2):
	"""
		Computes arrival time of direct, and two critically refracted waves from three layer models

		x: offset (array)
		v1: velocity of 1st layer (float)
		v2: velocity of 2nd layer (float)
		v3: velocity of 3rd layer (float)
		z1: thickness of 1st layer (float)
		z2: thickness of 2nd layer (float)

	"""
	direct = 1./v1*x
	theta1 = np.arcsin(v1/v2)    
	theta2 = np.arcsin(v2/v3)
	ti1 = 2*z1*np.cos(theta1)/v1
	ti2 = 2*z2*np.cos(theta2)/v2
	xc1 = 2*z1*np.tan(theta1)
	xc2 = 2*z2*np.tan(theta2)
	act1 = x > xc1
	act2 = x > xc2   
	ref1 = np.ones_like(x)*np.nan
	ref1[act1] = 1./v2*x[act1] + ti1    
	ref2 = np.ones_like(x)*np.nan
	ref2[act2] = 1./v3*x[act2] + ti2 + ti1
	t0 = 2.*z1/v1
	refl1 = np.sqrt(t0**2+x**2/v1**2)

	return np.c_[direct, ref1, ref2, refl1]

def evalExactφ(sqp,sqpnext,sqpdesired):

	a=sqp/sqpdesired
	b=sqpnext/sqpdesired
	if(abs(a+1) > float("1e-12") and (a*a + b*b - 1)>0):

		φ1=2*np.arctan(     ( b -  (a*a + b*b - 1)**0.5 )/(a+1)    )
		φ2=2*np.arctan(     ( b +  (a*a + b*b - 1)**0.5 )/(a+1)    )

		if(np.isnan(φ1)):
			print("Ok nan found, locating butter chicken")
			print( (a*a + b*b -1) )
			print( (b + (a*a + b*b -1)**0.5 ) )

		φ1=np.arcsin(np.sin(φ1))
		φ2=np.arcsin(np.sin(φ2))
		if(abs(φ1)<abs(φ2)):
			print(φ1,φ2)
			return φ1
		else:
			print(φ2,φ1)
			return φ2
	else:
		print("Glaba, something went globular :P")
		return 0

    def reconstruct_cluster_angle(self, events, index_group):
        """Reconstruct angles from a single event"""

        c = 3.00e+8

        t = [int(events[u]['ext_timestamp']) for u in index_group]
        stations = [events[u]['station_id'] for u in index_group]

        dt1 = t[0] - t[1]
        dt2 = t[0] - t[2]

        r1, phi1 = self.cluster.calc_r_and_phi_for_stations(stations[0], stations[1])
        r2, phi2 = self.cluster.calc_r_and_phi_for_stations(stations[0], stations[2])

        phi = arctan2((dt2 * r1 * cos(phi1) - dt1 * r2 * cos(phi2)),
                      (dt2 * r1 * sin(phi1) - dt1 * r2 * sin(phi2)) * -1)
        theta1 = arcsin(c * dt1 * 1e-9 / (r1 * cos(phi - phi1)))
        theta2 = arcsin(c * dt2 * 1e-9 / (r2 * cos(phi - phi2)))

        e1 = sqrt(self.rel_theta1_errorsq(theta1, phi, phi1, phi2, r1, r2))
        e2 = sqrt(self.rel_theta2_errorsq(theta2, phi, phi1, phi2, r1, r2))

        theta_wgt = (1 / e1 * theta1 + 1 / e2 * theta2) / (1 / e1 + 1 / e2)

        if theta_wgt < 0:
            theta_wgt *= -1
            phi += pi
            phi = (phi + pi) % (2 * pi) - pi

        return theta_wgt, phi

 def calc_total_reflection(self,alpha,n1,n2):
     beta=90.-(np.arcsin(n1/n2*np.sin(alpha)))*180/np.pi
     thetaC=np.arcsin(n1/n2)*180/np.pi
     if beta>thetaC:
         return True
     else:
         return False

def sphericalLawOfSines(angle1, side1, angle2, side2=None):
 '''Law of sines. 4 possible inputs, angle1 and side1 must be dms arrays and 
    either angle2 or side2 must also be dms arrays, with the other as None
    This function will determing the missing side/angle from the other three
    parameters
 '''
 #if side2 must be found
 if side2 == None:
  #do math, notice the conversion to/from radians
  angle1 = dms2deg(angle1)*np.pi/180
  angle2 = dms2deg(angle2)*np.pi/180
  side1  = dms2deg( side2)*np.pi/180
  side2  = np.arcsin(np.sin(angle2)*np.sin(side1)/np.sin(angle1))*180/np.pi

  return deg2dms(side2)

 #if angle2 must be found
 if angle2 == None:
  #do math, notice the conversion to/from radians
  angle1 = dms2deg(angle1)*np.pi/180
  side1  = dms2deg( side1)*np.pi/180
  side2  = dms2deg( side2)*np.pi/180
  angle2 = np.arcsin(np.sin(side2)*np.sin(angle1)/np.sin(side1))*180/np.pi

  return deg2dms(angle2)

 #tell the user that they fucked up
 else:
  print('Error, law of sines overconstrained, returning 0')
  return np.array([0,0,0])

    def test_arcsin_no_warning_on_unscaled_quantity(self):
        a = 15 * u.kpc
        b = 27 * u.pc

        with warnings.catch_warnings():
            warnings.filterwarnings('error')
            np.arcsin(b/a)

    def _transform_parametersraw(self, pars, in_format, out_format):
        """Convert parameter values from in_format to out_format.

           Also restores parameters to a preferred range if it permits multiple
           values that correspond to the same physical result.

           Parameters
           pars: Sequence of parameters
           in_format: A format defined for this class
           out_format: A format defined for this class

           Defined Formats
           internal: [position, parameterized width-squared, area]
           pwa: [position, full width at half maximum, area]
           mu_sigma_area: [mu, sigma, area]
        """
        temp = np.array(pars)

        # Do I need to change anything?  The internal parameters may need to be
        # placed into the preferred range, even though their interpretation does
        # not change.
        if in_format == out_format and in_format != "internal":
            return pars

        # Convert to intermediate format "internal"
        if in_format == "internal":
            # put the parameter for width in the "physical" quadrant [-pi/2,pi/2],
            # where .5*(sin(p)+1) covers fwhm = [0, maxwidth]
            n = np.floor((temp[1]+np.pi/2)/np.pi)
            if np.mod(n, 2) == 0:
                temp[1] = temp[1] - np.pi*n
            else:
                temp[1] = np.pi*n - temp[1]
            temp[2] = np.abs(temp[2]) # map negative area to equivalent positive one
        elif in_format == "pwa":
            if temp[1] > self.maxwidth:
                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(temp[1], self.maxwidth)
                raise SrMiseTransformationError(emsg)
            temp[1] = np.arcsin(2.*temp[1]**2/self.maxwidth**2-1.)
        elif in_format == "mu_sigma_area":
            fwhm = temp[1]*self.sigma2fwhm
            if fwhm > self.maxwidth:
                emsg = "Width %s (FWHM) greater than maximum allowed width %s" %(fwhm, self.maxwidth)
                raise SrMiseTransformationError(emsg)
            temp[1] = np.arcsin(2.*fwhm**2/self.maxwidth**2-1.)
        else:
            raise ValueError("Argument 'in_format' must be one of %s." \
                              % self.parformats)

        # Convert to specified output format from "internal" format.
        if out_format == "internal":
            pass
        elif out_format == "pwa":
            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)
        elif out_format == "mu_sigma_area":
            temp[1] = np.sqrt(.5*(np.sin(temp[1])+1.)*self.maxwidth**2)/self.sigma2fwhm
        else:
            raise ValueError("Argument 'out_format' must be one of %s." \
                              % self.parformats)
        return temp