Well, I could take a guess at this one. The elevation angle for any object is given by: Sin(E) = Cos(D) Cos(L) - Sin(D) Sin(L) Sin(H) where E = elevation angle, D = sun's declination, L = observers latitude and H = the hour angle at the time of the observation. For sunrise and sunset, E = 0.0 by definition, so for a given Latitude on the earth, there will be two values for H that will solve this equation: sunrise and sunset. Now, the sun's declination changes continuously throughout the year and depends on the day (season). A very rough formula for this is: D = 24.0 sin(d) where d = 360*( t )/365 degrees, and t is the number of days since the Spring Equinox. This gives declinations of 0.0, +24.0, 0.0 and -24.0 at t = 0.0, 90, 180 and 270 days for the Sun at the Spring, Summer, Autumn and Winter ( Northern Hemisphere) ecliptic positions which are at RIght Ascension 12h, 6h, 0h and 18h. OK, we now know D, and you have to provide L for each of the latitudes that you will be using to plot this curve, that leaves H.
H depends on the Right Ascension of the observer and the local time which we take to be Universal Time. This is related to the Longitude of the observer on the Earth at sunrise and sunset. RA = 12 - 24*t/365 where t is the number of days since the Spring Equinox and remember that when RA goes negative you have to add 24.0 to it. As a guess, I would say that H = RA - Long with RA converted into degrees by multiplying by 360/24. and with Longitude running westward from 0 - 360 degrees beginning at the prime meridian in Greenwich.
Now, you should be able to use these formulae like this. First select a day of the year ( season) and compute D and RA. Then set E = 0.0. You now have an equation that looks like: A*Cos(L) = B*Sin(L) *Sin(H). On your Mercator map, you then find the pair of latitudes and longitudes ( sunrise and sunset) that solve the equation.
I haven't had any time to test any of this, so good luck!