The goal of this research was to develop an analytic low-thrust trajectory design method. The term low thrust stems form the fact that current electric propulsion systems can only produce low amounts of thrust. Nevertheless, these systems can deliver very high amounts of ?V using little propellant, which makes them very suitable for interplanetary flight. The orbital mechanics regarding low-thrust interplanetary flight were considered and it was found that for first-order design it is acceptable to neglect all perturbing forces with respect to the Sun's gravitational attraction and the (low) thrust force. In addition to that, a look was taken at the characteristics of electric propulsion and feasible thrust (acceleration) levels were obtained. Since only low thrust levels can be achieved by electric propulsion, long periods of thrusting are required to obtain sufficient change in velocity. This continuous thrusting feature of electrically-propelled spacecraft highly complicates the dynamics of such spacecraft. As a consequence, the design and optimization of low-thrust trajectories is very difficult, and therefore preliminary design of such trajectories is very important. Such first-order designs can give a fast overview of all feasible trajectories and can, moreover, be used as an initial guess for more refined trajectory optimization. For generating and optimizing first-order low-thrust trajectories, analytical methods are very suited, since they are very fast and exact. Several analytical low-thrust methods have been developed recently and most of them are shape-based; the trajectory is assumed to have a specific shape which satisfies the boundary conditions and subsequently the thrust profile can be obtained from the dynamics. These shape-based methods are popular for preliminary design, because of their simplicity and ability to obtain feasible trajectories. The currently known shape-based methods were investigated to obtain know-how about their functioning and to find their strengths and weaknesses. It was noticed that regarding the solving of boundary conditions improvements could be made and that most shaping methods lack flexibility. Inspired by the use of velocity hodographs for the computation of non-perturbed transfer orbits, a novel low-thrust trajectory design method was developed based on shaping the velocity of spacecraft during the transfer. For the shaping of the velocity components, velocity functions were used which consist of a sum of simple base functions. These base functions can be integrated analytically, such that the change in position can be obtained analytically. Since the velocity is shaped, the conditions on initial and final velocity can be solved very easily and exactly. In addition, the boundary conditions on position can be solved exactly without the need of iterative computations. Next to the parameters required to satisfy the boundary conditions, extra parameters can be added to make the design and optimization of trajectories more flexible. Two different methods have been developed; one which shapes the velocity as a function of time and another one which shapes as a function of the polar angle. For the computation of the required ?V , and for the computation of the polar angle in the time-driven method, numerical integration is required. Four simple integration methods were tested and an RK4 integrator was found to perform best based on accuracy and computation speed. Since the derivative of the ?V and the polar angle (the thrust acceleration and angular velocity, respectively) can be computed as a function of time or polar angle, their values can be computed prior to the integration. As a result, RK4 integration simplifies to Simpson's rule and the integration can be done much faster. Both hodographic-shaping methods have been verified by numerically propagating the initial spacecraft state vector using thrust profiles found by the shaping methods. Both the time-driven and polar-angle-driven method were found to function well, and the obtained trajectories and thrust profiles agreed almost exactly. In order to obtain minimum-?V trajectories, the free parameters in the velocity functions were optimized. Both a Nelder-Mead (NM) and Differential Evolution (DE) optimizer were tested for this purpose. NM was found to be equally robust, but much faster than DE, and therefore NM was used. In addition, initial guesses for the free parameters were used to speed up the optimization. Finally, the search for the optimal departure date and TOF was done by stepping through the flight window using a grid. Both hodographic-shaping methods were tested for missions to Mars, the near-Earth asteroid 1989ML, comet Tempel 1 and Mercury. For missions to Mars and 1989ML the two methods gave comparable results. The time-driven method found better trajectories to Tempel 1 (a far outer target), whereas the polar-angle-driven method found better ones to Mercury (an inner target). In comparison with the spherical and pseudo-equinoctial shaping methods and DITAN (the benchmark), the hodographic-shaping methods performed well, especially for mission to Mars and to asteroid 1989ML. The near-optimal results came at the cost of computational speed. The best results were found using 6 free parameters requiring require a computation time of 2 s per trajectory on average. The lowest-order solutions (0 DoF), on the other hand, required only 1.6 ms per trajectory on average. Finally, the lowest-order solutions were found to be unable to obtain near-optimal trajectories, however, they were able to indicate the correct regions of low ?V in the flight windows and to give a good indication of the required ?V and thrust acceleration.