Thank you, Mikko for doing such a deep dive into my program!
From the manual I read the UliSpeed cover "5-DOF", or 5 degrees of freedom, only yaw missing from a full 6-DOF VPP. Most VPPs are still only 2-DOF (I believe the ORC is), boatspeed & heel, ignoring even the effect of leeway. You have the rudder angle as an input, and also estimate the canoe body side force/drag, so would it not be possible to add the 6th DOF, by iterating the rudder angle and underwater Mz against the sail Mz, in a similar fashion as you have leeway as an independent parameter within every iteration?
Theoretically, it would be possible. I use XFOIL to determine the aerodynamic coefficients of the local profile of the sail. The lift coefficient is well predicted by XFOIL, the drag is still estimated closely enough for a VPP, but the position of the center of effort can be erroneous, especially if the flow is separated. Calculating the balance would add a huge computational load for a questionable result. On the other hand, I am interested in the geometric optimization of the boat. If the first prototype on the water would require a rudder angle that is too large, one would change the position of the mast or of the keel, but would not substantially change the general geometry of the boat.
UliSpeed is also the only VPP I have seen that truly trims sails to maximize speed, with real world parameters such as sheeting angle, twist and camber. Camber only for the jib, though, are you defining the mainsail shape to maintain an as linear as possible downwash distribution for the sailplan?
To create the database for the mainsail with XFOIL requires the variation of the Reynolds-number, the sheeting angle, the mast diameter and the angle of attack, each with a spacing that is close enough to allow interpolation. Right now, the database comprises 1680 lift and 1680 drag coefficients and 1680 values for the camber. Each operating point uses the camber that gives the best lift to drag ratio for the input parameters. If the camber of the mainsail were added as a further input parameter, the database would be larger by the factor 8 and the runtime would increase similarly. The question, if it is possible to produce such a mainsail, is the weak point of my method.
From what I have seen, the proposed sail shapes perform very well, when analysed with CFD, even very different shapes providing quite similar results, as you state in the manual. On my wish list in the OUT-result file would be the mainsail and the jib CL & CD, easily integrated from the given Cl & Cd distributions. This would highlight the role of the main vs. the jib, as the jib CD is often close to zero, very much different from the ORC sailcoefficients, for instance.
I will do that and add the forces created by the individual sails. However, this information could be disturbing to the novice (not for you

). The jib operates in a flow field that is distorted by the main and vice versa. If each sail were used alone, the forces would be quite different. Max Munk described this in his biplane theory in 1923.
Aviation and Aerodynamics History: Max M Munk https://aerodyn.org/munk/ (stagger theorem).
In my program, the induced velocities caused by the free and bound vortices of the sails and of their mirror images are calculated using the Biot-Savart-law. In addition, the effect of the streamline curvature is taken into account. As a result, the flow field for each sail is quite different from the undisturbed onset flow. The resultant force vector at the jib is rotated forward, whereas the vector at the main is rotated backwards. One consequence is for example: If you leave everything unchanged and just reef the mainsail, the force on the jib will change and its drag-component will increase.
Thanks again, Uli