訂正
× for (n = 0; n < 5 && fabs(ckrad) >= EPSILON * sinrad; n++) {
○ for (n = 0; n < 500 && fabs(ckrad) >= EPSILON * sinrad; n++) {

これでもradが大きくなるとずれてくるのがわかるはず