def METHOD_NAME(self,zl_min=0,zl_max=1100,n_zl=10,log_zl=False,pk_func=None,
            zs1=[1100],p_zs1=[1],zs2=[1100],p_zs2=[1],
            pk_params=None,cosmo_h=None,l=np.arange(3001)):
    if not cosmo_h:
        cosmo_h=self.cosmo_h
    if log_zl:#bins for z_lens. 
        zl=np.logspace(np.log10(max(zl_min,1.e-4)),np.log10(zl_max),n_zl)
    else:
        zl=np.linspace(zl_min,zl_max,n_zl)
    ds1=cosmo_h.comoving_transverse_distance(zs1)
    ds2=cosmo_h.comoving_transverse_distance(zs2) 
    dl=cosmo_h.comoving_transverse_distance(zl)
    prefac=1.5*(cosmo_h.H0/c)**2*cosmo_h.Om0
    sigma_c1=prefac*(1+zl)*dl*(1-np.multiply.outer(1./ds1,dl))
    sigma_c2=prefac*(1+zl)*dl*(1-np.multiply.outer(1./ds2,dl))
    clz=self.cl_z(z=zl,l=l,cosmo_h=cosmo_h,pk_params=pk_params,pk_func=pk_func)
    dzl=np.gradient(zl)
    dzs1=np.gradient(zs1) if len(zs1)>1 else 1
    dzs2=np.gradient(zs2) if len(zs2)>1 else 1
    cl_zs_12=np.einsum('ji,ki,il',sigma_c2,sigma_c1*dzl,clz)#integrate over zl..
    cl=np.dot(p_zs2*dzs2,np.dot(p_zs1*dzs1,cl_zs_12))
    cl/=np.sum(p_zs2*dzs2)*np.sum(p_zs1*dzs1)
    f=l*(l+1.)/(l+0.5)**2 #correction from Kilbinger+ 2017
    cl*=f**2
    return l,cl