Was doing some testing for high latitude and the date where it has 2 sunset is off by one when compare to timeanddate.com and GML sunrise sunset date.

so after doing some digging to the website algorithm that code uses.

The solar transit approximation in J is using slightly off coefficient.

0.0069 should be 0.0068 for Earth which correspond to J2 in equation 35.

After updating the coefficient, the calculation for double sunset date is matching what is shown in timeanddate.com.