ได้ลองดูซอสโค้ดโปรแกรมภาษาซีของ Swiss Ephemeris หรือยังครับ

ในตัวสูตรที่ใช้กับ Quick Basic ถ้าตัดเฉพาะคำนวณดาวทิพย์จะเป็นทำนองนี้ครับ ยังไม่รู้ว่าโพสต์แล้วมันจะรักษาช่องเว้นวรรคข้างหน้าให้หรือเปล่านะครับ

ค่าคงที่ของดาวทิพย์

DATA CU,5526,82.4854457,6.400766885e-3,-2.755224344e-5,
 DATA 33,62,85,100,106,101,87,65,41,16,-6,-20,-25,-19,-4,19,49,81
 DATA HA,-1965,61.23270944,-2.772376054e-03,9.452816256e-6
 DATA 0,26,50,69,80,94,79,68,50,30,10,-8,-19,-24,-20,-8,10,34,60
 DATA ZE,5794,48.3720491,-4.285060477e-3,3.864569431e-5
 DATA 0,22,43,58,68,71,67,56,41,24,6,-9,-19,-23,-20,-10,6,25,47
 DATA KR,669,42.16960545,-3.117970292e-3,2.221915822e-5
 DATA 0,20,38,53,61,64,60,51,37,21,5,-9,-18,-22,-20,-11,3,21,42
 DATA AP,7931,37.64981827,-2.907258758e-3,2.009076057e-5
 DATA 0,20,36,49,57,59,54,47,32,18,2,-11,-19,-22,-20,-11,2,17,38
 DATA AD,-864,35.28683374,-5.003003022e-3,2.73466561e-5
 DATA 0,17,33,46,53,55,52,42,31,16,1,-11,-19,-22,-20,-11,2,17,35
 DATA VU,3039,32.74576048,-6.107844531e-4,1.531796472e-6
 DATA 0,18,32,44,52,53,50,42,29,15,1,-10,-17,-21,-18,-11,-1,14,33
 DATA PO,9660,29.23595125,-3.149299571e-3,1.789787019e-5
 DATA 0,17,31,43,48,50,47,40,29,15,2,-8,-18,-20,-17,-10,1,14,28

ส่วนการคำนวณ

DAOTIP:
 READ x$, A, b, c, d
 FOR I = 0 TO 18
  READ A(I)
 NEXT I
   z = Y0 - 1890: jj = TOJJ(A, b, c, d, z)
 IF jj >= 270 AND LS > jj THEN GOTO DaotipOut
 IF (m0 > 3 OR LS < 30) AND LS > jj THEN GOTO DaotipOut
 IF (m0 > 30 OR LS < 30) AND jj >= 270 THEN GOTO DaotipOut
   z = z - 1: jj = TOJJ(A, b, c, d, z)
DaotipOut:
 j = LS - jj: j = TON(j): CALL TOXX(A(), j, XX)
 j = jj + XX: j = TON(j): GOTO Sompusdao
Sompusdao:
 j = INT(j * 60 + .5): L(K) = j / 60
 P = L(K): Mi = j - INT(P) * 60: dm = INT(P) + Mi / 100
 IF K <= 14 THEN PRINT "       "; x$; : PRINT USING " ###.##"; dm; : x = L(K): CALL SOMPUS(x)
 IF K > 14 THEN LOCATE (K - 8), 45: PRINT x$; : PRINT USING " ###.##"; dm; : x = L(K): CALL SOMPUS(x)

ฟังก์ชันที่ใช้

SUB SOMPUS (x)
 CALL SDM(x, S, d, m)
 PRINT USING "   ##  ##  ## "; S; d; m
END SUB

FUNCTION TOJJ (A, b, c, d, z)
 j = A + b * z + c * z ^ 2 + d * z ^ 3 + .5
 j = INT(j) / 60: TOJJ = TON(j)
END FUNCTION

FUNCTION TON (x)
 TON = x - INT(x / 360) * 360
END FUNCTION

SUB TOXX (A(), j, XX)
 x = j
 SELECT CASE x
  CASE IS <= 120: N = 0: m = 6: GOTO Cal
  CASE IS <= 240: N = 6: m = 12: GOTO Cal
  CASE IS > 240: N = 12: m = 18
 END SELECT
Cal:
 y = 0
 FOR j% = N TO m
 P = 1
  FOR I% = N TO m
   IF I% = j% THEN GOTO FIN
   P = P * (x - I% * 20) / (j% * 20 - I% * 20)
FIN:
  NEXT I%
  y = y + P * A(j%)
 NEXT j%
 XX = INT(y) / 60
END SUB

ขอบคุณอาจารย์อีกครั้งครับ

จำได้ว่าในกระทู้เก่าๆ เคยมีคนขอดูซ้อสโค้ดรุ่นนี้ตัวเต็ม ตอนนั้นผมปฏิเสธไปว่า ในซ้อสโค้ดมันมีชื่อผู้เขียนโปรแกรมอยู่ ซึ่งที่จริงมันเป็นการดีที่ให้เครดิตเขา แต่เกรงจะเป็นเหมือนเผยแพร่โดยไม่ได้รับอนุญาต คงต้องขอดูเหตุผลประกอบอื่นๆ ด้วยครับ

ฟังก์ชันที่ Error นี่คงเป็นการแปลงค่าทศนิยมเป็นราศี-องศา-ลิบดา ถ้ายังไงพอจะเขียนเองเพิ่มเติมไปก่อนได้ไหมครับ