This paper introduces a novel methodology for robust intersatellite ranging and attitude measurement planning via a case study of two low Earth orbit small satellites flying in formation. The relative position model express the curvilinear coordinates dynamics of the Follower in the Leader centered local frame, and includes deterministic effects of the J2 acceleration. The relative attitude model express the linearized dynamics of the Euler angles from the Leader centered local frame to the Follower body frame and includes the Gravity-Gradient torque. The atmospheric air density is modeled as a stochastic process perturbing both accelerations and torques. Its means and variances are calibrated via Monte-Carlo simulations of a high-fidelity Truth model. Measured intersatellite ranges and relative attitudes are processed along a finite period of time via a Kalman filter to produce estimates of a twelve states vector of relative position and attitude and of their rates. The estimation error covariance matrix at the final time possess a known upper bound function of the observability and noise controllability Gramians. First an optimal ranging strategy is developed by finding the ranging noise variance profile that minimizes the upper bound. The minimization is subject to a time integral constraint on the noise variance that is derived from a laser energy budget. This results in a remarkably short sequence of ranging epoch times. Then that problem is extended to the maximization of the upper bound with respect to the atmospheric density variance profile. The density variance is assumed to satisfy a time integral constraint, derived from finite energy considerations and quantified through Monte-Carlo simulations. The problems are solved through numerical iterative schemes and their solutions consist of very few ranging acquisition times, atmospheric density impulses, along with the corresponding optimized intensities. The combined solution to both problems thus provides a robust laser ranging planning over a given window of time. Robustness is here in terms of guaranteed performances on the estimation error covariance matrix at the final time: its upper bound is designed for the worst-case atmospheric density profile along the scheduled trajectory. The proposed design has got appealing practical features such as: very few ranging epoch times compared with continuous ranging, for a given overall laser energy budget; very few process noise impulses, compared with continuous process noise, for a given overall noise energy. Extensive simulations are performed that verify the algorithms convergence and validate the proposed approach. Extensive Monte-Carlo simulations are also presented to compare the performances of Kalman filters designed via the proposed approach with others based on continuous measurement and process noise.