I need to draw maps on a computer screen.

From Wikipedia, I found that for a given latitude, the "y" value for a Mercator projection is given by...

y = ln(tan(Latitude) + (1 / Cos(Latitude)))

... so I wrote a program to calculate "y" for latitudes from 0.0 to 89.9.

Now from a "sit back and look at it" position, the "y" value increases the further you move from the equator, so if I have a 1 degree square cell, when I plot it, with say 100 pixels for 1 degree of longitude, the scaling due to the projection should mean that I multiply that same 100 in longitude by the factor I get from the formula above, and that I should end up with a figure larger then 100, i.e. that "y" calculated above should be greater then 1 for all values of latitude.

It is not. I have to set latitude = 50.0 before the result of the above formula passes 1. I am multiplying the latitude in degrees by Pi/180.0 to convert the degrees to Radians. I am using 3.14159265358979323846 as Pi. All of my variables are either declared to be double precision floating points, or explicitly stated with a decoimal point to force float operation by the compiler.

Is the formula wrong? Am I doing something silly?