! 異種外射方位図法による海岸線（地球儀）の表示プログラム
! 実行時に、海岸線データを記入した外部ファイル "coast.txt" を要求。
! ダウンロード後、ソースプログラムと同じディレクトリに置いてください。
! 外部ファイル内、１行ごとに「経度,緯度」。線を切るところは「0,0」
! 海岸線データは、以下より取得したもの。
! http://rimmer.ngdc.noaa.gov/coast/  WCL DESIGNED FOR 1:5,000,000
!
! （十進BASIC ソース）：　全文コピーして実行してください。
!
!                            2010.02.
!                            2015.04. 改訂 
!                            S. Nakajima
!   十進BASIC は、以下のURIより入手してください。
!   http://hp.vector.co.jp/authors/VA008683/


! ３次元の極座標を３次元デカルト座標に変換する定義式
! ただし、R, T, F は r, θ, Φ
! R は地球の半径
! T は北極で 0、赤道は PI/2、南極で PI
! F は経度（西経は0〜-PI、東経は0〜+PI）
DEF X(R,T,F) = R * SIN(T) * COS(F)
DEF Y(R,T,F) = R * SIN(T) * SIN(F)
DEF Z(R,T,F) = R * COS(T)

! ３次元の x,y,z 座標を、二次元上で表示する際の表示座標の定義式。
LET TH = PI/2*0.5      ! z軸を手前に倒す角度
LET FI = PI/2*0.28     ! z軸周りの回転
DEF XX(X,Y,Z) =  -X *SIN(FI)          +Y *COS(FI)          
DEF YY(X,Y,Z) =  -X *COS(FI)*SIN(TH)  -Y *SIN(FI)*SIN(TH)   +Z *COS(TH)

! 画面準備
ASK bitmap SIZE h,v
LET a = h/v
SET WINDOW -11*a,11*a,-11,11

! 半径一定 → 球の半径を定義。
LET R = 10


DO   
   CLEAR   ! 画面の消去
    
   ! 背景の表示
   SET AXIS COLOR 15
   DRAW GRID 
    
    
   INPUT PROMPT "自転軸を手前に倒す【南半球 -1 〜 1 北半球】" : TT
   INPUT PROMPT "自転方向に回す 【グリニッジ基準: 0 〜 1 /周】" : FF   
   LET TH = PI/2 *TT     ! z軸を手前に倒す角度
   LET FI = PI*2 *(-FF)  ! z軸周りの回転
    
   ! 地球儀部分の背景を塗りつぶす
   SET AREA COLOR 0   ! 白
   DRAW DISK WITH SCALE(r)
    
   ! 影線の表示
   SET LINE COLOR 15
   FOR s = 0 TO 0.99 STEP 0.01
      DRAW circle WITH SCALE(r*SQR(1-s^2))
   NEXT s   
    
   SET LINE COLOR 8   ! 緯経線の色。8 は灰色。    
   ! 輪郭線の表示 
   DRAW circle WITH SCALE(r)
    
   ! 緯度線の表示
   FOR T = 0 TO PI STEP PI/30
      FOR F = FI TO FI+2.1*PI STEP PI/30
      ! 手前側を判定し、それ以外は描画しない
         IF SIN(-T)*COS(F-FI)*COS(TH) - COS(-T)*SIN(TH) <0 THEN   
            PLOT LINES : XX(X(R,T,F),Y(R,T,F),Z(R,T,F)), YY(X(R,T,F),Y(R,T,F),Z(R,T,F));
         ELSE
            PLOT LINES
         END if
      NEXT F
      PLOT LINES
   NEXT T  
    
   ! 経度線の表示
   FOR F = 0 TO PI*2.1 STEP PI/10
      FOR T = PI/60 TO 59/60*PI STEP PI/120
         IF SIN(-T)*COS(F-FI)*COS(TH) - COS(-T)*SIN(TH) <= 0 THEN  
            PLOT LINES : XX(X(R,T,F),Y(R,T,F),Z(R,T,F)), YY(X(R,T,F),Y(R,T,F),Z(R,T,F));
         END if
      NEXT T
      PLOT LINES
   NEXT F
    
   ! 赤道の表示
   SET LINE COLOR 4 ! 赤
   LET T = PI/2 
   FOR F = 0 TO 2.1*PI STEP PI/120
      IF SIN(-T)*COS(F-FI)*COS(TH) - COS(-T)*SIN(TH) < 1/100 THEN   
         PLOT LINES : XX(X(R,T,F),Y(R,T,F),Z(R,T,F)), YY(X(R,T,F),Y(R,T,F),Z(R,T,F));
      ELSE
         PLOT LINES
      END if
   NEXT F
   PLOT LINES
    
    
   ! 海岸線データの読み込みと表示
   ! 外部ファイル。１行に「経度,緯度」。線を切るところは「0,0」
    
   SET LINE COLOR 1  ! 黒
    
   OPEN #1 : NAME "coast.txt" ,ACCESS INPUT
   DO
      INPUT #1 , IF MISSING THEN EXIT DO : A,B
      IF A=0 THEN
         PLOT LINES
      ELSE
         LET T = (90-B) *PI/180
         LET F = A *PI/180
         IF SIN(-T)*COS(F-FI)*COS(TH) - COS(-T)*SIN(TH) <0 THEN   
            PLOT LINES : XX(X(R,T,F),Y(R,T,F),Z(R,T,F)), YY(X(R,T,F),Y(R,T,F),Z(R,T,F));
         ELSE
            PLOT LINES
         END IF
      END IF 
   LOOP
    
   CLOSE #1
   PLOT LINES
    
   LET DT = 0   ! 再描画までの待機時間（秒）
   LET DELAY$ = ""
   IF DT > 0 THEN LET DELAY$ = STR$(DT) & "秒経過後に"
   PAUSE "OK を押すと、" & DELAY$ & "画面をクリアし、再描画します。" 
   WAIT DELAY DT
    
LOOP
END