Atmospheric correction of Earth Observation data is one of the most critical steps in the data processing chain of a satellite mission for successful remote sensing applications. Atmospheric Radiative Transfer Models (RTM) inversion methods are typically preferred due to their high accuracy. However, the execution of RTMs on a pixel-per-pixel basis is impractical due to their high computation time, thus large multi-dimensional look-up tables (LUTs) are precomputed for their later interpolation. To further reduce the RTM computation burden and the error in LUT interpolation, we have developed a method to automatically select the minimum and optimal set of nodes to be included in a LUT. We present the gradient-based automatic LUT generator algorithm (GALGA) which relies on the notion of an acquisition function that incorporates (a) the Jacobian evaluation of an RTM, and (b) information about the multivariate distribution of the current nodes. We illustrate the capabilities of GALGA in the automatic construction and optimization of MODerate resolution atmospheric TRANsmission (MODTRAN) LUTs for several input dimensions. Our results indicate that, when compared to a pseudo-random homogeneous distribution of the LUT nodes, GALGA reduces (1) the LUT size by $\sim$75\% and (2) the maximum interpolation relative errors by 0.5\% It is concluded that automatic LUT design might benefit from the methodology proposed in GALGA to reduce computation time and interpolation errors.