This is the code that finally works.

thetashould be the angle between optic axis and the point on the image, so tan(theta) will be r/f assuming ris the distance of the 2D point (on the image) from the center of the image. Looks like I got the core concept wrong in the original post

img = cv2.imread(<impath>)[:,:,::-1]
H, W, _ = img.shape
cX, cY = W//2, H//2 #7, 5 
f = min(cX, cY)
mapX = np.zeros((H, W))
mapY = np.zeros((H, W))
for x in range(W):
    for y in range(H):
        dx, dy = (x - cX), (y - cY)
        r = np.sqrt(dx**2 + dy**2)
        phi = np.arctan2(dy, dx)
        theta = np.arctan(r/f)
        rnew = f*theta
        xnew = cX + rnew*np.cos(phi)
        ynew = cY + rnew*np.sin(phi)
        mapX[H - 1 - int(ynew), int(xnew)] = H - 1 - y
        mapY[H - 1 - int(ynew), int(xnew)] = x
distorted = cv2.remap(np.array(img, "uint8"), np.array(mapY, "float32"), np.array(mapX, "float32"), interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT)