-- main.lua
-- ****** Comment *****************************************************
-- Revision: 1.0
-- Issued  : Feb.15, 2026
-- Auther  : Tanzai
-- Contents: This script controlls the model, "H3_Spacecraft" on X-Plane 12.
--           This Plugin shall be stored in /H3_Spacecraft_X.X/plugins/ folder under airplane
--           folder.
-- Notice 1: The package may be redistributed for non commercial purposes
--           only; derivatives MAY NOT be redistributed.
-- Notice 2: Do not resize X-Plane window during flight to prevent incorrrect
--           display. At any rate you can restore the window by pressing a menu 
--           "Reload theCurrent Airplane (Skip Art Reload)" under "Developer".
-- *******************************************************************
-- Custom slider setting
--   showPart[10] = Boosters smoke/jet display    0:hide  1:show
--   showPart[11] = Stage 1 jet display           0:hide  1:show
--   showPart[12] = Stage 2 jet display           0:hide  1:show
--   showPart[13] = Re-entry jet display          0:hide  1:show

setAircraftPanelRendering(false)
set3DRendering(false)
setInteractivity(false)

-- Variables setting  ---------------------------------------------------
defineProperty("startIsOn", 0)
defineProperty("simStage", 0)
defineProperty("resetIsOn", 0)
defineProperty("mapIsOn", 0)              -- 1: draws map
defineProperty("isInternal", 0)
defineProperty("t3", 0)
defineProperty("tDisp", 0)
defineProperty("Note1", " ")              -- 飛行中画面の説明書き
defineProperty("crrVelo", 0.0)            -- 現在速度(m/s)
defineProperty("curLon", 0.0)             -- 現在経度(度)
defineProperty("curLat", 0.0)             -- 現在緯度(度)
defineProperty("kmH", 0.0)                -- 表示高度(km)
defineProperty("portId", "")              -- 目的地空港ID
defineProperty("portName", "NOT FOUND")   -- 目的地空港ID
defineProperty("isKeyReady", 0 )          -- 目的地空港ID
defineProperty("lonD", 0 )                -- 目的地空港経度
defineProperty("latD", 0 )                -- 目的地空港緯度
defineProperty("isStartReady", 0)         -- Start押下可能
defineProperty("isInLimit", 0)            -- Cabin内に視点がある
defineProperty("driftX", 0.0)            -- Cabin内のX方向浮動量
defineProperty("driftY", 0.0)            -- Cabin内のY方向浮動量
defineProperty("driftZ", 0.0)            -- Cabin内のZ方向浮動量

addSearchResourcesPath ( moduleDirectory .. "/configuration/" )
hideBooster = createGlobalPropertyi("H3/hide/Booster", 0)    -- 1:hide
hideChuteL = createGlobalPropertyi("H3/hide/ChuteL", 0)      -- 1:hide
hideChuteS = createGlobalPropertyi("H3/hide/ChuteS", 0)      -- 1:hide
hideFirstS  = createGlobalPropertyi("H3/hide/FirstStage", 0) -- 1:hide
hideFairing = createGlobalPropertyi("H3/hide/Fairing", 0)    -- 1:hide
hideSecondS = createGlobalPropertyi("H3/hide/SecondStage",0) -- 1:hide
hideShip1 = createGlobalPropertyi("H3/hide/Ship1",0)         -- 1:hide
hideShip1B = createGlobalPropertyi("H3/hide/Ship1B",1)       -- 1:hide
hideShip2 = createGlobalPropertyi("H3/hide/Ship2",1)         -- 1:hide
points1= createGlobalPropertyfa("H3/root/points1" , 73 ) -- points(x1,y1,..)
points2= createGlobalPropertyfa("H3/root/points2" , 73 ) -- points(x1,y1,..)
LauncherObj = sasl.loadObject( "H3_Launch_pad.obj" )
H3_coundown = sasl.al.loadSample ( "H3_countdown.wav" )
H3_sound = sasl.al.loadSample ( "H3_sound.wav" )
as = globalPropertyf("sim/flightmodel/position/true_airspeed") -- m/s 速度
elv = globalPropertyd("sim/flightmodel/position/elevation")    -- n, m 高度
lat = globalPropertyd("sim/flightmodel/position/latitude") -- n, deg 緯度
lon = globalPropertyd("sim/flightmodel/position/longitude")-- n, deg 経度
quat = globalPropertyfa("sim/flightmodel/position/q")      -- y, 四元数
vx = globalPropertyf("sim/flightmodel/position/local_vx")  -- y, m/s 東が正
vy = globalPropertyf("sim/flightmodel/position/local_vy")  -- y, m/s 真上が正
vz = globalPropertyf("sim/flightmodel/position/local_vz")  -- y, m/s 南が正
planeX = globalPropertyd("sim/flightmodel/position/local_x") -- 東を正
planeY = globalPropertyd("sim/flightmodel/position/local_y") -- 真上が正
planeZ = globalPropertyd("sim/flightmodel/position/local_z") -- 南が正
vRoll = globalPropertyf("sim/flightmodel/position/P") -- ロ-ル回転速度(deg/s)
vPitch= globalPropertyf("sim/flightmodel/position/Q") -- ピッチ回転速度(deg/s)
vYaw  = globalPropertyf("sim/flightmodel/position/R") -- ヨー回転速度(deg/s)
headX = globalPropertyf("sim/graphics/view/pilots_head_x")
headY = globalPropertyf("sim/graphics/view/pilots_head_y")
headZ = globalPropertyf("sim/graphics/view/pilots_head_z")
headThe = globalPropertyf("sim/graphics/view/pilots_head_the") -- Pilots pitch
headPsi = globalPropertyf("sim/graphics/view/pilots_head_psi") -- Pilots heading
sim_time = globalPropertyf("sim/operation/misc/frame_rate_period")
wheelZoom = globalPropertyi("sim/graphics/view/mouse_wheel_zoom_external")
isExternal = globalPropertyi("sim/graphics/view/view_is_external")
Booster_release = createGlobalPropertyf("H3/Booster/release", 0)
First_release   = createGlobalPropertyf("H3/First/release", 0)
Second_release  = createGlobalPropertyf("H3/Second/release", 0)
Fairing_release = createGlobalPropertyf("H3/Fairing/release", 0)
showPart = globalPropertyia("sim/cockpit2/switches/custom_slider_on")
animTime = sasl.createTimer()
inner_view = sasl.findCommand ( "sim/view/3d_cockpit_cmnd_look" )
viewSet = sasl.findCommand ( "sim/view/circle" )
-- zoomOnExt = globalPropertyi("sim/graphics/view/mouse_wheel_zoom_external")
oldVelo = {0.0,0.0,0.0}
xFactor = {}
yFactor = {}
zFactor = {}
fuelMass = {}
iniMass = {}
totalMass = {}
points1 = {}
points2 = {}
pResult = {}

-- Subwindow setting  ------------------------------------------
-- winW = globalPropertyi("sim/graphics/view/window_width")
winH = globalPropertyi("sim/graphics/view/window_height")

local Map = contextWindow {
--	position 		= { 0, 0, get(winW), get(winH)-80},
	position 		= { 0, 0, 1198,509},
	noResize 		= true,
	visible 		= true,
	noBackground 	= true,
	layer			= SASL_CW_LAYER_FLOATING_WINDOWS,
	noDecore 		= true,
	noMove 			= true,
--	components 		= {mapArea{position={0,0,get(winW),get(winH)-80},},},
	components 		= {mapArea{position={0,0, 1198,509},},},
}

local infoArea = contextWindow {
    position		= { 0, get(winH)-80, 1250, 80},
    noResize		= false,
    visible			= true,
	proportional	= false,
    noBackground    = true,
	layer			= SASL_CW_LAYER_FLOATING_WINDOWS,
    noDecore        = true,
    noMove          = true,
    components		= {simDraw {position={0, 0,1250, 75}, bg = {1,1,1,1}},
 					   simInput{position={0, 0,1200, 34},},},
}

-- Function to set rocket parameters  -----------------------------
function setRocket()
	targetH = 120                            -- 目標高度 120 km とする
	-- Spec   Booster   1st    2nd  Orbiter
	tStart  ={   0  ,  110  ,  185 ,  0   }  -- 燃焼開始時間(整数) sec
	firingT ={ 105  ,   70  ,  300 ,  0   }  -- 燃焼時間(整数) sec

	-- 最初の6秒間は迎角90°となります。
	-- Boost終了からFirst点火まで5秒以上開けること。
	-- Firstの燃焼時間は12秒以上の事。
	-- First終了からSecond点火まで5秒以上開けること。
	-- Second燃焼中に高度がtargetHを超えたらに迎角を0°に変更する。
	-- Second燃焼中に速度がtargetVを超えたら停止する。
	-- 燃焼時間(firingT)の上限設定は、firingT={ 105, 153, 370, 0 }

	-- Spec    Booster   1st     2nd  Orbiter
	gLsp    ={ 2279.3, 4165.0, 4390.4,  0.0 }  -- 有効排気速度(=g*Lsp) m/s
	dryMass ={   37.4,   12.5,    3.0,  3.0 }  -- 空筐重量 ton
	fuelRate={  1.082, 0.7092, 0.03108, 0.0 }  -- 燃料消費速度 ton/sec

	for m = 1,4 do
		fuelMass[m] = fuelRate[m] * firingT[m]
		iniMass[m] = dryMass[m] + fuelMass[m]
	end
	totalMass[1] = iniMass[1] + iniMass[2] + iniMass[3] + iniMass[4]
	totalMass[2] = iniMass[2] + iniMass[3] + iniMass[4]
	totalMass[3] = iniMass[3] + iniMass[4]
end

-- Other functions  ----------------------------------------------
function setOrient(psi,theta,phi)
	psSR = math.sin(psi/2)        -- Heading N is 0 (rad)
	psCR = math.cos(psi/2)
	thSR = math.sin(theta/2)      -- pitch XZ plane is 0 (rad)
	thCR = math.cos(theta/2)
	phSR = math.sin(phi/2)        -- roll  (rad)
	phCR = math.cos(phi/2)
	q0 =  psCR*thCR*phCR + psSR*thSR*phSR
	q1 =  psCR*thCR*phSR - psSR*thSR*phCR
	q2 =  psCR*thSR*phCR + psSR*thCR*phSR
	q3 = -psCR*thSR*phSR + psSR*thCR*phCR
	set(quat, {q0,q1,q2,q3})
end

function setSpeed(n)                 -- Setting speed by Tsiolkovsky's equation
	xFactor[n] = -math.sin(pitchThe) * math.sin(headPsiR)
	yFactor[n] = math.cos(pitchThe)
	zFactor[n] = math.sin(pitchThe) * math.cos(headPsiR)
	-------------
	massRatio= 1.0 - fuelMass[n]*(get(t3)-tStart[n])/firingT[n]/totalMass[n]
--
--   speedDrag = cd * a * psl * exp(-beta*h) * V^2 /2
--	cd = 1.2                                       抵抗係数 (ジェミニは1.55)
--  a = 21.237                                     機体正面面積(m2) (ジェミニは4.143)
--	Psl = 1.225                                    海面空気密度(kg/m3)
--	Beta = 1/6.7                                   スケール高度の逆数(1/km)
--	speedDrag = 15.6 * exp(-h/6.7) * V^2 
--
	deltaVelo= gLsp[n]*math.log(1/massRatio)
--	deltaVelo= gLsp[n]*math.log(1/massRatio) - 15.6/math.exp(get(elv)/6700)*oldVelo[n]^2
	vInX = (oldVelo[n] + deltaVelo) * xFactor[n]         -- 東西方向の速度
	vInY = (oldVelo[n] + deltaVelo) * yFactor[n] - 9.8   -- 垂直方向係数、重力を加味する
	vInZ = (oldVelo[n] + deltaVelo) * zFactor[n]         -- 南北方向の速度
	set(vx, vInX)
	set(vy, vInY)
	set(vz, vInZ)
	vDelta = (vInX^2 + vInY^2 + vInZ^2)^0.5
	set(crrVelo, vDelta)                                 -- crrVelo: 速度m/s
end

function setFreeFall(n)                                  -- Setting free fall speed
	vInY = vInY - 9.8                                    -- 重力を加味
	set(vy, vInY)
	vDelta = (vInX^2 + vInY^2 + vInZ^2)^0.5
	set(crrVelo, vDelta )                -- crrVelo: 1秒毎の速度
end

function stayRocketThere()
	set(vx, 0.0)
	set(vy, 0.0)
	set(vz, 0.0)
end

function stopRotation()
	set(vRoll, 0.0)
	set(vPitch, 0.0)
	set(vYaw, 0.0)
end

function trimEndD(valD)
	if valD > 180 then valD = valD - 360 end
	if valD < -180 then valD = valD + 360 end
	return valD
end

function trimEndR(valR)
	if valR > math.pi then valR = valR - 2*math.pi end
	if valR < -math.pi then valR = valR + 2*math.pi end
	return valR
end

function getAirportInfo( newId )
	NavAidID = sasl.findNavAid(nil,newId,nil,nil,nil,NAV_AIRPORT )
	type,la,lo,destH,freq,heading,id,name,isInDSF = sasl.getNavAidInfo(NavAidID)
	if newId == id and newId ~= "" then
		set(portName, name)
		set(lonD, lo)
		set(latD, la)
	else                                   -- set Mt.Fuji as default
		set(portName, "Top of Mt.Fuji")
		set(portId, "????")
		id = "Mt.Fuji"
		set(lonD, 138.730487)
		set(latD, 35.362473)
		lo = 138.730487
		la = 35.362473
		destH = 3750
	end
	return id, name, la, lo, destH
end

function getOrbParam2(lonAD, latAD, lonBD, latBD)
-- 地球中心と地点A,Bを含む軌道面のパラメータを求める。
	RA1, incD = getOrbPlane(lonAD, latAD, lonBD, latBD)
	RA = RA1
   
	if RA ~= 0 then                       -- 経度をRAだけ回転し昇交点を仮のX軸に合わせる。
		lonAD2 = trimEndD(lonAD - RA)
		lonBD2 = trimEndD(lonBD - RA)
		RA1, incD = getOrbPlane(lonAD2, latAD, lonBD2, latBD)
		if incD > 0 then               -- 昇交点が仮のX軸と一致しているので、この座標を採用する。
			orbAD = getOrbD(math.rad(lonAD2), math.rad(latAD))
			orbBD = getOrbD(math.rad(lonBD2), math.rad(latBD))
			nBar,deltaBarR = getNBar(lonAD2,latAD,lonBD2,latBD,orbAD,orbBD)
		elseif incD < 0 then -- 昇交点が仮のX軸のマイナス側と一致しているので、座標を更に180度回転する。
			RA = trimEndD(RA + 180)
			lonAD3 = trimEndD(lonAD2 + 180)
			lonBD3 = trimEndD(lonBD2 + 180)
			RA1, incD = getOrbPlane(lonAD3, latAD, lonBD3, latBD)
			orbAD = getOrbD(math.rad(lonAD3), math.rad(latAD))
			orbBD = getOrbD(math.rad(lonBD3), math.rad(latBD))
			nBar,deltaBarR = getNBar(lonAD3,latAD,lonBD3,latBD,orbAD,orbBD)
		else                            -- incD=0 即ち軌道面が赤道面に一致した場合
			RA = 0
			orbAD = lonAD
			orbBD = lonBD
		end
	else                             -- RA=0 の場合
		RA = 0
		incD = 0
		orbAD = getOrbD(math.rad(lonAD), math.rad(latAD))
		orbBD = getOrbD(math.rad(lonBD), math.rad(latBD))
	end
	print("Final RA=",RA,"incD=",incD,lonAD,latAD,"orbAD=",orbAD,"orbBD=",orbBD)
	return RA, incD, orbAD, orbBD, nBar, deltaBarR
end

function getOrbPlane(lonAD, latAD, lonBD, latBD)
-- 地点A,B を通る軌道面の軌道傾斜角(incD)、昇交点赤径(RA)を求める。
-- 軌道面方程式  ax + by + cz = 0
	Ax = math.cos(math.rad(latAD))*math.cos(math.rad(lonAD))
	Ay = math.cos(math.rad(latAD))*math.sin(math.rad(lonAD))
	Az = math.sin(math.rad(latAD))
	Bx = math.cos(math.rad(latBD))*math.cos(math.rad(lonBD))
	By = math.cos(math.rad(latBD))*math.sin(math.rad(lonBD))
	Bz = math.sin(math.rad(latBD))
	a = Ay*Bz - Az*By
	b = Az*Bx - Ax*Bz
	c = Ax*By - Ay*Bx
	RA1 = trimEndD(math.deg( math.atan( -a/b )))
	incD = math.deg(math.atan( -b/c ))
	return RA1, incD
end

function getOrbD(lonR, latR)
-- 昇交点赤経にX軸を合わせた座標で地点(lonR,latR)の軌道面内角(orbD)を求める。
	orbD = math.deg(math.acos( math.cos(latR)*math.cos(lonR)))
	if latR < 0 then orbD = -orbD end
	orbD = trimEndD(orbD) 
	return orbD
end

function getOrbR(lonR, latR)
-- 昇交点赤経にX軸を合わせた座標で地点(lonR,latR)の軌道面内角(orbR)を求める。
	orbR = math.acos( math.cos(latR)*math.cos(lonR))
	if latR < 0 then orbR = -orbR end
	orbR = trimEndR(orbR) 
	return orbR
end

function getLonLat2(orbR, incR, RA)
-- 昇交点赤経(RA)、軌道傾斜角(incR)、面内角(orbR)からその地点の緯度、経度、を求める。
	newLat = math.deg(math.asin(math.sin(orbR)*math.sin(incR)))
	newLon = math.deg(math.atan(math.tan(orbR)*math.cos(incR))) + RA
	if orbR >  math.pi/2 then newLon = newLon + 180 end
	if orbR < -math.pi/2 then newLon = newLon - 180 end
	newLon = trimEndD( newLon )
	return newLon, newLat
end

function getNBar(lonA, latA, lonB, latB, orbA, orbB )
-- 2地点A,Bの経度、緯度、面内角から軌道描画の分解数(nBar)や分解角(deltaBarR)を求める。
	if orbA*orbB > 0 then
		if orbA < orbB then
			deltaOrb = math.abs(orbB - orbA)
		else
			deltaOrb = 360 - math.abs(orbA + orbB)
		end
	elseif orbA > 0 then
		deltaOrb = 360 + orbB - orbA
	elseif lonA < 0 then
		deltaOrb = orbB - orbA
	end
	stepAng = 10
	nBar = math.floor(deltaOrb/stepAng)
	deltaBarR = math.rad(deltaOrb/nBar)
	return nBar, deltaBarR
end

function getRootMap2(lonAD,latAD,lonBD,latBD,RAAN,incD,orbAD,orbBD,nBar,deltaBarR)
	orbR = math.rad(orbAD)
	iChange = 0
	
	for i = 1, 73, 2 do
		if i == 1 then
			points1[i] = lonAD
			points1[i+1] = latAD
		elseif i < 1+2*nBar then
			orbR = trimEndR(orbR + deltaBarR)
			newLon,newLat = getLonLat2(orbR, math.rad(incD), RAAN)
			points1[i]   = newLon
			points1[i+1] = newLat
			headPsiR = getPsi(newLon, newLat, RAAN)
			if iChange==0 and points1[i]<points1[i-2] then
				iChange = i-2                          -- iChange番目までは最初の曲線
			end
		elseif i >= 1+2*nBar then
			points1[i] = lonBD
			points1[i+1] = latBD
		end
	end

	if iChange == 0 then
		for j = 1, 73, 2 do
			points1[j]   = 599 + points1[j]*598/180            -- 経度α
			points1[j+1] = 212 + points1[j+1]*304/90           -- 緯度β
			points2[j]   = -180
			points2[j+1] = -90
		end
	else
		for j = 1, 73, 2 do
			if j <= iChange then
				points2[j]   = points1[iChange+2]
				points2[j+1] = points1[iChange+3]
			elseif j <= 2*nBar+1 then
				points2[j]   = points1[j]
				points2[j+1] = points1[j+1]
				points1[j]   = points1[iChange]
				points1[j+1] = points1[iChange+1]
			elseif j > 2*nBar+1 then
				points1[j]   = points1[iChange]
				points1[j+1] = points1[iChange+1]
				points2[j]   = points2[2*nBar+1]
				points2[j+1] = points2[2*nBar+2]
			end
		end
			
		for j = 1, 73, 2 do
			points1[j]   = 599 + points1[j]*598/180            -- 経度α
			points1[j+1] = 212 + points1[j+1]*304/90           -- 緯度β
			points2[j]   = 599 + points2[j]*598/180            -- 経度α
			points2[j+1] = 212 + points2[j+1]*304/90           -- 緯度β
		end
	end
end

function drawLauncher(ps)
	iniX = get(planeX)
	iniY = get(planeY) - 1
	iniZ = get(planeZ)
	Launcher = sasl.createInstance(LauncherObj)
	sasl.setInstancePosition(Launcher,iniX,iniY,iniZ,0,ps,0)
end

function getPsi(lonA, latA, RA)
-- 昇交点赤経(RA)、ある地点の緯度(lonA)、経度(latA)からその地点の飛行方位(psi)を求める。
	lonR = trimEndR(math.rad(lonA - RA))
	latR = math.rad(latA)
	xR = math.cos(latR)*math.cos(lonR)
	yR = math.cos(latR)*math.sin(lonR)
	zR = math.sin(latR)
	Psi = math.acos(math.abs(xR*zR)/((yR^2 + xR^2)^0.5)/((zR^2 + yR^2)^0.5))
	k0 = trimEndD(lonA - RA - 90)
	if k0 > 0 and k0 < 180 then
		Psi = math.pi - Psi
	end
	return Psi
end	

-- sim_time後の進行角vrad2、速度vero2、高度ht2、進入角ang2を求めるFunction
-- dV/dt = -0.00232834 *exp(-h/6.7) * V^2 - 396675070/r~2 * sinγ
-- dγ/dt= (V/r - 396675070/r~2/V) * cosγ
-- h = h + V * sinγ / 1000

function getRadH(x1, x2, h, dt)
	-- x1:速度(m/s), x2:進入角(rad), h:高度(km), dt:時間間隔(秒)
    k11 = dt * func1(x1, x2, h)
    k12 = dt * func2(x1, x2, h)
    k21 = dt * func1(x1+0.5*k11, x2+0.5*k12, h)
    k22 = dt * func2(x1+0.5*k11, x2+0.5*k12, h)
    k31 = dt * func1(x1+0.5*k21, x2+0.5*k22, h)
    k32 = dt * func2(x1+0.5*k21, x2+0.5*k22, h)
    k41 = dt * func1(x1+k31, x2+k32, h)
    k42 = dt * func2(x1+k31, x2+k32, h)
    x1 = x1 + (k11/6 + k21/3 + k31/3 + k41/6)
    x2 = x2 + (k12/6 + k22/3 + k32/3 + k42/6)
    
    dRad = x1*math.cos(x2)*dt/6360000        -- dt秒後の移動角(rad)
    h2 = h + x1*math.sin(x2)/1000            -- dt秒後の高度(km)
	gVal = -func1(x1, x2, h) * dt /9.8       -- dt秒後の重力加速度(G)
	return dRad, h2, gVal, x1, x2
end

function gValue(h)                   -- get G value at height h(km)
    gV = 396675070/(h + 6360)^2
	return gV
end

function func1(x1, x2, h)          -- get dV/dt (accelaration)
    dVdT = -0.00232834*math.exp(-h/6.7)*x1^2 - gValue(h)*math.sin(x2)
	return dVdT
end
 
function func2(x1, x2, h)          -- get dGamma/dt
    dGdT = (x1/(h+6360)/1000 - gValue(h)/x1) * math.cos(x2)
	return dGdT
end

function iniValSet()              -- Function to set initial value
	t2 = -6
	pitchThe = 0                  -- 直立時を0radとする迎角(モデルが水平の時は-1.57079となる)
	set(Note1, "")
	set(showPart, 0, 10)          -- No show Boosters smoke/jet
	set(showPart, 0, 11)          -- No show First Stage jet
	set(showPart, 0, 12)          -- No show Second Stage jet
	set(showPart, 0, 13)          -- No show Backfire
	Action2 = 0
	Action4 = 0
	Action5 = 0
	Action8 = 0
	set(Booster_release, 0)
	set(First_release, 0)
	set(Second_release, 0)
	set(Fairing_release, 0)
	set(hideBooster,0)
	set(hideFirstS, 0)
	set(hideFairing,0)
	set(hideSecondS, 0)
	set(hideShip1, 0)
	set(hideShip1B, 1)
	set(hideShip2, 1)
	set(hideChuteL, 1)
	set(hideChuteS, 1)
	state = 0
	state1 = 0
	state2 = 0
	count1 = 0                 -- for print out counter
	count2 = 0                 -- for setting at reaching targetspeed
	count3 = 0
	count4 = 0
	count5 = 0
	nLine = 24                 -- マップの軌道線分の数
	set(isKeyReady, 0)
	set(isStartReady, 0)
	preKey = 1
	vRad = 0.0
	endVelo = 7832.03          -- 周回速度 m/s
	set(portId, "")
	set(headPsi, 0.0)
	set(headThe, 0.0)
	set(driftX, 0.0)
	set(driftY, 0.0)
	set(driftZ, 0.0)
	set(isInternal, 0)
	prePsi = 0.0               -- set initial view psi
	preThe = 0.0               -- set initial view the
	tDesc = 0
	preV = 1
end

-- Main control ----------------------------------------------------
set(simStage, 0)
set(startIsOn, 0)
set(resetIsOn, 0)
iniValSet()                                     -- Initialize values
set(mapIsOn, 1)
setRocket()                                     -- Set rocket parameters
sasl.setCSMode(1)
set(wheelZoom, 1)
--[[
-- script for Ship2 debug
	set(mapIsOn, 0)
	prePsi = 0.0            -- set initial view psi
	preThe = 0.0            -- set initial view the
	pitchThe = 0.0
	sasl.setCSMode(1)
	set(hideFirstS, 1)
	set(hideShip2, 0)
	sasl.commandOnce( inner_view )
]]
function update()
	-----------  Setting before launch  -----------
	if get(simStage) == 0 then                  -- 軌道設定画面
		stayRocketThere()
		stopRotation()

		if preKey == 1 and get(isKeyReady) == 0 then
			set(curLon, get(lon))
			set(curLat, get(lat))
			id,name,la,lo,destH = getAirportInfo(get(portId)) -- Airport確認
			print("Flight to:", id, name, la, lo, destH)
			set(isStartReady, 1)
			RAAN,incD,orbAD,orbBD,nBar,deltaBarR = getOrbParam2(get(curLon),get(curLat),lo,la)
			print(get(curLon),get(curLat),lo,la,RAAN,incD,orbAD,orbBD,nBar,deltaBarR)
			getRootMap2(get(curLon),get(curLat),lo,la,RAAN,incD,orbAD,orbBD,nBar,deltaBarR)
			set(Note1, "by         min.")
			
			-- 軌道上にある軌道角 angTop を求める。軌道到達までの軌道角は13.39377°、
			-- 突入後の軌道角は25.17726°、合計38.55149°
			if orbBD<0 and orbAD>0 then orbBD = orbBD + 360 end
			angTop = math.rad(orbBD - orbAD - 38.55149)  -- 軌道上にある期間の軌道角
			if angTop < 0 then
				angTop = angTop + 2*math.pi
			end
			set(tDisp, math.floor((angTop/0.0012052069 + 905 )/60))
			orbToDesR = math.rad(orbBD - 25.17726)       -- 突入開始角
			if orbToDesR < -math.pi then orbToDesR = orbToDesR + 2*math.pi end
			print(orbAD,math.deg(angTop),math.deg(orbToDesR),orbBD,tDisp)
			
			preKey = 0
		end
		
		if get(isKeyReady)==1 then              -- 目的地の入力が可能になった。
			preKey = 1                          -- 目的地の入力が現在可能
		end
		
		if get(startIsOn) == 1 then             -- "Liftoff" pressed 
			iniLon = get(curLon)
			iniLat = get(curLat)
			iniElv = get(elv)                            -- (m)
			destLonR = math.rad(lo - RAAN)
			destLatR = math.rad(la)
			destOrbR = getOrbR(destLonR, destLatR)
			headPsiR = getPsi(iniLon, iniLat, RAAN)
			iniPsi = math.deg(headPsiR)
			pitchThe = 0.0
			rollPhi = 0.0
			sasl.startTimer( animTime )
			---------
			set(isStartReady, 0)
			set(simStage, 1)
			set(startIsOn, 0)
			iniValSet()                         -- Initialize values
			set(mapIsOn, 0)
			sasl.commandOnce( viewSet )
			setRocket()                         -- Set rocket parameters
		end
		
		if get(resetIsOn) == 1 then             -- "Reset" pressed 
			set(startIsOn, 0)
			set(resetIsOn, 0)
			sasl.placeUserAtLocation(iniLat,iniLon,iniElv,iniPsi,0.0)
			iniValSet()                         -- Initialize values
			sasl.commandOnce( viewSet )
			setRocket()                         -- Set rocket parameters
		end
	end
	
	-----------  Control during launching  -----------
	if get(simStage) == 1 then                      -- Simulation starts
		t1 = sasl.getElapsedSeconds( animTime )
		if math.floor(t1)-7 >= t2 then              -- 1秒毎に以下を実施
			t2 = t2 + 1
			set(t3, t2)
			set(curLon, get(lon))                   -- 経度 deg
			set(curLat, get(lat))                   -- 緯度 deg
			set(kmH, get(elv)/1000)

			if t2 == -5 then
				sasl.al.playSample( H3_coundown )
			elseif t2 < 0 then
				stayRocketThere()
				set(Note1, "Countdown starts")
			elseif t2 == 0 then                     -- Booster fires
				tEnd1 = tStart[1] + firingT[1]
				oldVelo[1] = 0.0
				setSpeed(1)
				set(vx, 0.0)
				set(vz, 0.0)
				sasl.al.playSample( H3_sound )
				set(showPart, 1, 10)                -- Boosters jet on
				set(Note1, "Lifting up by Boosters until T= 105")
				curLonR = math.rad(get(curLon) - RAAN)
				orbRad = getOrbR(curLonR, math.rad(get(curLat)))
				print("Liftoff at", t2, math.deg(orbRad))
			elseif t2 < tEnd1 then
				setSpeed(1)
			elseif t2 == tEnd1 then                -- Booster stops
				setSpeed(1)
				set(showPart, 0, 10)               -- Boosters jet off
				Action2 = 1          -- Booster release animation starts 5s
			elseif t2 < tStart[2] then
				setFreeFall(1)
				set(Note1, "Boosters released")
 			elseif t2 == tStart[2] then            -- 1st Stage starts
				tEnd2 = tStart[2] + firingT[2]
				oldVelo[2] = get(crrVelo)
				setSpeed(2)
				set(showPart, 1, 11)               -- Stage 1 jet on
				Action4 = 1           -- Animation of Fairing release 20s
				set(Note1, "1st Stage engine lifting up until T= 180")
			elseif t2 < tEnd2 then
				setSpeed(2)
			elseif t2 == tEnd2 then                -- 1st Stage ends
				setSpeed(2)
				set(showPart, 0, 11)               -- Stage 1 jet off
				Action5 = 1                    -- animation Stage 1 release
				set(Note1, "1st Stage engine cutoff")
			elseif t2 < tStart[3] then
				setFreeFall(2)
			elseif t2 == tStart[3] then            -- 2nd Stage starts
				tEnd3 = tStart[3] + firingT[3]
				oldVelo[3] = get(crrVelo)
				setSpeed(3)
				set(showPart, 1, 12)               -- Stage 2 jet on
				set(Note1, "2nd Stage engine running")
			elseif t2 < tEnd3 then
				setSpeed(3)
				if get(kmH) > targetH then     -- Reach target height
					if count5 == 0 then
						set(Note1, "Keeping height, Increasing speed until 7830 m/s")
						count5 = 1
					end
					if get(crrVelo) >= endVelo then
						set(simStage, 2)
						sasl.setCSMode(1)
						Action8 = 1                    -- Release 2nd Stage
						set(Note1, "Orbiting around the Earth until T=")
						set(showPart, 0, 12)           -- Stage 2 jet off
						curLonR = math.rad(get(curLon) - RAAN)
						orbRad = getOrbR(curLonR, math.rad(get(curLat)))
						vRad1 = get(crrVelo)/(get(kmH)+6378.14)/1000  -- 角速度(rad/s)
						tOrb = angTop/vRad1             -- 軌道上にある時間 sec
						set(tDisp, get(t3) + math.floor(tOrb))
						print("Orbiting starts from", t2, math.deg(orbRad), vRad1)
					end
				end
			elseif t2 >= tEnd3 then      -- 2nd Stage stops out of orbit
				set(Note1, "Failed to reach target height and speed")
			end	
		end
	
		if Action2 == 1 then             -- Remove Boosters
			state = get(Booster_release) + get(sim_time)/5
			set(Booster_release, state)
			if state >= 1 then
				set(hideBooster,1)
				Action2 = 0 
				state = 0
			end
		end

		if Action4 == 1 then             -- Animation of Fairing release
			state1 = get(Fairing_release) + get(sim_time)/5
			set(Fairing_release, state1)
			if state1 >= 1 then
 				set(hideFairing,1)
				Action4 = 0
				state1 = 0
			end
		end
		
		if Action5 == 1 then              -- animation of Stage 1
			state = get(First_release) + get(sim_time)/10
			set(First_release, state)
			if state >= 0.5 then
				set(hideFirstS,  1)       -- 1st Stage removed
				Action5 = 0
				state = 0
			end
		end
		pitchThe = -1.55334*get(elv)/120000
		setOrient(headPsiR,pitchThe,rollPhi)
		stopRotation()
	end

	---------  Control during orbiting  ---------------- 
	if get(simStage) == 2 then                        -- Orbiting phase
		t1 = sasl.getElapsedSeconds( animTime )
		if math.floor(t1)-7 >= t2 then                -- 秒毎に以下を実施
			t2 = t2 + 1
			set(t3, t2)
			set(curLon, get(lon))                       -- 経度 deg
			set(curLat, get(lat))                       -- 緯度 deg
			phiR = getPsi(get(curLon), get(curLat), RAAN)
			set(vx, get(crrVelo)*math.sin(phiR))        -- 東西方向の速度
			set(vy, 0.0)                                -- 垂直方向の速度
			set(vz, -get(crrVelo)*math.cos(phiR))       -- 南北方向の速度
			rollPhi = phiR - headPsiR
		end
		
		orbRad = orbRad + vRad1*get(sim_time)                   -- 180°で補正しない軌道角
		if orbRad > math.pi then
			orbRad = orbRad - 2*math.pi
		end
		newLon,newLat = getLonLat2(orbRad, math.rad(incD), RAAN)
		newX,newY,newZ = sasl.worldToLocal(newLat,newLon,get(kmH)*1000)
		set( planeX, newX )
		set( planeY, newY )
		set( planeZ, newZ )
		
		if orbRad > orbToDesR then                    -- 予定時間経過したので突入する
			orbRad = orbToDesR + 0.0006
			set(simStage, 3)
			set(Note1, "Descending into Earth to touch at T=")
			set(tDisp, get(t3) + 655 )
			set(isExternal,1)
			sasl.commandOnce( viewSet )
			sasl.commandBegin( viewSet )
			sasl.commandEnd( inner_view )
			pitchThe = 1.570795
			x1 = get(crrVelo)   -- 初期速度(m/s)
			x2 = -0.0279252     -- 初期経路角 γ = -1.6°
			dt = 1              -- 計算間隔(1 sec)
			dRad = x1*math.cos(x2)*dt/6360000    -- 突入中の面内角速度(rad)
			mtH = get(elv)
			print(t2,x1,x2,get(kmH),"Descend from",math.deg(orbRad),"to",orbBD)
		end
		
		if Action8 == 1 then                        -- Release 2nd stage
			state = state + get(sim_time)/4
			set(Second_release, state)
			if state >= 1 then
				set(hideSecondS, 1)
				set(hideShip1, 1)
				set(hideShip2, 0)
				Action8 = 0
				state = 0
			end
		end
		setOrient(headPsiR,pitchThe,rollPhi)
		stopRotation()
	end

	---------  Control during re-enry  ---------------- 
	if get(simStage) == 3 then                        -- re-enry phase
		-- 次のorbRadを求め、その緯度経度に機体を移動する。
		orbRad = orbRad + dRad*get(sim_time)           -- 現在のorbRad
		if orbRad > math.pi then orbRad = orbRad - 2*math.pi end
		newLon,newLat = getLonLat2(orbRad, math.rad(incD), RAAN)
		mtH = mtH + x1*math.sin(x2)*get(sim_time)   -- 現在高度(m)
		newX,newY,newZ = sasl.worldToLocal(newLat,newLon,mtH)
		set( planeX, newX )
		set( planeY, newY )
		set( planeZ, newZ )
		if mtH < destH + 3000 then         -- 目的地上空3000mなので着地体勢に入る
			set(hideChuteS, 0)             -- パラシュート小を開く
			set(simStage, 4)
			set(Note1, "Descending by drogue parachute")
			print(t2,"above 3000 m", mtH, math.deg(orbRad))
		end
			
		t1 = sasl.getElapsedSeconds( animTime )
		if math.floor( t1 - 6 ) == t2 then              -- 秒毎に以下を実施
			t2 = t2 + 1
			tDesc = tDesc +1 
			set(t3, t2)
			
			-- 現在の緯度経度速度高度進入角からdt(1秒)後のo速度高度進入角を求める
			dRad, h2, gVal, x1, x2 = getRadH(x1, x2, get(kmH), dt)
			set(kmH, h2)
		
			-- フレーム毎に次の地点の情報設定
			set(curLon, get(lon))                       -- 経度 deg
			set(curLat, get(lat))                       -- 緯度 deg
			set(crrVelo, math.ceil(10*x1)/10)
						
			-- 現在位置での速度方向指定
			phiR = getPsi(get(curLon), get(curLat), RAAN)
			set(vx,  x1*math.cos(x2)*math.sin(phiR)) -- 東西方向速度
			set(vy,  x1*math.sin(x2))                -- 垂直方向速度
			set(vz, -x1*math.cos(x2)*math.cos(phiR)) -- 南北方向速度
			pitchThe = 1.570795 + x2
			rollPhi  = -phiR+headPsiR
			
			if t2%10 == 0 then             -- 10秒毎に以下を出力
				angG = math.deg(x2)
				angO = math.deg(orbRad)
				print(t2,get(crrVelo),angG,get(kmH),angO,get(curLon),get(curLat))
			end
			
			if tDesc == 180 then set(showPart,1,13)      -- backfire on
			elseif tDesc == 280 then
				set(hideShip2, 1)
				set(hideShip1B, 0)
			elseif tDesc == 370 then set(showPart,0,13)  -- backfire off	
			end
		end
		setOrient(headPsiR,pitchThe,rollPhi)
		stopRotation()
	end
	
		---------  Control at touchdown  ---------------- 
	if get(simStage) == 4 then                       -- touchdown phase
		-- フレーム毎にorbRadを求め、その緯度経度に機体を移動する。
		if mtH > destH + 200 then
			mtH = mtH - 40.0*get(sim_time)           -- 高度200までは40m/sで降下
		else
			mtH = mtH - 4.0*get(sim_time)            -- 高度200以下は4m/sで降下
			if count4 == 0 then
				set(hideChuteS, 1)                   -- パラシュート小を閉じる
				set(hideChuteL, 0)                   -- パラシュート大を開く
				sasl.commandBegin( viewSet )
				set(Note1, "Descending by main parachutes")
				count4 = 1
			end
		end
		newX,newY,newZ = sasl.worldToLocal(get(latD),get(lonD),mtH)
		set( planeX, newX )
		set( planeY, newY )
		set( planeZ, newZ )
		
		if mtH <= destH +24 then                     -- 機体が着地した
			set(vy, 0.0)                             -- 垂直方向速度
			set(hideChuteL, 1)                       -- パラシュート大を閉じる開く
			set(simStage, 5)
			set(Note1, "Landed on terrain at T=")
			set(tDisp, get(t3))
			print(t2,"Spacecraft hits terrain", mtH,get(lon),get(lat))
		end
				
		t1 = sasl.getElapsedSeconds( animTime )
		if math.floor( t1 - 6 ) == t2 then              -- 1秒毎に以下を実施
			t2 = t2 + 1
			set(t3, t2)
			
			-- 現在の緯度経度速度高度進入角からdt秒後の速度高度進入角を求める
			set(curLon, get(lonD))                       -- 目的地経度 deg
			set(curLat, get(latD))                       -- 目的地緯度 deg
			if mtH > destH + 200 then set(crrVelo, -40.0)
			else
				set(crrVelo, -4.0)
			end
			set(kmH, mtH/1000)
			if t2%10 == 0 then
				angG = -90                               -- 進入角(°)
				print(t2,get(crrVelo),angG,get(kmH),get(lonD),get(latD))
			end
			
			-- 現在位置での速度方向指定
			set(vx,  0.0)                            -- 東西方向速度
			set(vy, get(crrVelo))                    -- 垂直方向速度
			set(vz,  0.0)                            -- 南北方向速度
		end
		setOrient(0,0,0)                   -- (headPsiR,pitchThe,rollPh)
		stopRotation()
	end
	
		---------  Stay on ground  ---------------- 
	if get(simStage) == 5 then                       -- touchdown phase
		stayRocketThere()
		setOrient(0,0,0)                   -- (headPsiR,pitchThe,rollPh)
		stopRotation()
	end

	-- ----------- Cabin limit and Zero gravity effect  ----------
	if get(isExternal) == 0 then                      -- View is Internal
		set(isInLimit, 0)
		if (-20.8<get(headY) and get(headY)<-19.4) then
			if get(headZ) < 0.7 then
				wallR = get(headX)^2 + get(headZ)^2
				limitR = (1.55 + 0.25*(-19.4-get(headY)))^2
				if wallR < limitR then
					set(isInLimit, 1)
				end
			end
		end
	
		dragD = sasl.getCSDragDirection()
		dragV = sasl.getCSDragValue()
		if dragV ~= 0.0 then
			prePsi = prePsi + 0.005*dragV*math.sin(math.rad(90-dragD))
			preThe = preThe + 0.005*dragV*math.cos(math.rad(90-dragD))			
		end
		set(headPsi, prePsi)
		set(headThe, pitchThe + preThe)
		
		if sasl.getCSClickDown(MB_MIDDLE)~=1 then
			wheelC = sasl.getCSWheelClicks()
			if get(isInLimit) == 1 then                    -- View is in Cabin
				if wheelC == 0 then       -- Wheel not rotated
					set(headX, get(headX) + 0.02*get(driftX) )
					set(headY, get(headY) + 0.02*get(driftY) )
					set(headZ, get(headZ) + 0.02*get(driftZ) )
				end
				preV = 1
			end
			if get(isInLimit) ~= preV then      -- View is goint out of Cabin
				set(driftX, -get(driftX))
				set(driftY, -get(driftY))
				set(driftZ, -get(driftZ))
				preV = 0
				set(headX, get(headX) + get(driftX))
				set(headY, get(headY) + get(driftY))
				set(headZ, get(headZ) + get(driftZ))
			end
		else
			set(driftX, 0)
			set(driftY, 0)
			set(driftZ, 0)
		end
    end
end

-- default view setting  ------------------------------------------------
function onModuleInit ()
	drawLauncher( 0.0 )
end

function onAirportLoaded()
	sasl.commandBegin( viewSet )
	drawLauncher( 0.0 )
end