Given two $n \times n$ Hermitian matrices $Y$ and $Λ$, the Harish-Chandra-Itzykson-Zuber (HCIZ) distribution on the unitary group $\text{U}(n)$ is $e^{\text{tr}(UΛU^*Y)}dμ(U)$, where $μ$ is the Haar measure on $\text{U}(n)$. The density $e^{\text{tr}(UΛU^*Y)}$ is known as the HCIZ density. Random unitary matrices distributed according to the HCIZ density are important in various settings in physics and random matrix theory. However, the basic question of efficient sampling from the HCIZ distribution has remained open. We present two efficient algorithms to sample matrices from distributions that are close to the HCIZ distribution. The first algorithm outputs samples that are $ξ$-close in total variation distance and requires polynomially many arithmetic operations in $\log 1/ξ$ and the number of bits needed to encode $Y$ and $Λ$. The second algorithm comes with a stronger guarantee that the samples are $ξ$-close in infinity divergence, but the number of arithmetic operations depends polynomially on $1/ξ$, the number of bits needed to encode $Y$ and $Λ$, and the differences of the largest and the smallest eigenvalues of $Y$ and $Λ$. HCIZ densities can also be viewed as exponential densities on $\text{U}(n)$-orbits, and these densities have been studied in statistics, machine learning, and theoretical computer science. Thus our results have the following applications: 1) an efficient algorithm to sample from complex versions of matrix Langevin distributions studied in statistics, 2) an efficient algorithm to sample from continuous max-entropy distributions on unitary orbits, which implies an efficient algorithm to sample a pure quantum state from the entropy-maximizing ensemble representing a given density matrix, and 3) an efficient algorithm for differentially private rank-$k$ approximation, with improved utility bounds for $k>1$.