Preliminaries

BeginPackage["Interstellar`"];
(*  Interstellar package for making interstellar flight movies.  *)
(*  Copyright 1992 by Alfred M. Kaufman (AMK Development)  *)
(*  All rights reserved  *)
(*  Special thanks to Robert Managan of LLNL and John Stein of WRI *)

Contstants

c= 2.997925*10^10;  (*  LIGHT SPEED *)
g0 = 980.;  (*  Definition of 1 g *)
secsperyear = 3600.*24.*365.;
lightyear = c*secsperyear;
t0 = c/(g0*secsperyear);  (* Years  *)
begpt = .0004; (* light years *)
endpt = .0004; (* light years *)

Movie Frame

MovieFrame =
	{{GrayLevel[1.],Hue[.7],Rectangle[{-1.2,-1.2},{1.2,1.2}]},
	{GrayLevel[0.],Rectangle[{-1.,-1.},{1.,1.}]},
		{GrayLevel[1.],Rectangle[{-1.,1.03},{1.,1.17}]},
		{GrayLevel[1.],Rectangle[{-1.00,-1.17},{-.02,-1.03}]},
		{GrayLevel[1.],Rectangle[{.02,-1.17},{1.0,-1.03}]},
		{GrayLevel[1.],Rectangle[{-1.17,-1.0},{-1.03,1.0}]},
		{Hue[.7],Table[Line[{{-1.17,yy},{-1.03,yy}}],{yy,-1.0,1.0,.2}]},
		{GrayLevel[1.],Rectangle[{1.03,-1.0},{1.17,1.0}]},
		{Hue[.7],Table[Line[{{1.03,yy},{1.17,yy}}],{yy,-1.0,1.0,.2}]},
		{GrayLevel[1.],Rectangle[{-1.17,-1.17},{-1.03,-1.03}]},
		{GrayLevel[1.],Hue[.95],Text["v",{-1.1,-1.1}]},
		{GrayLevel[1.],Rectangle[{1.03,-1.17},{1.17,-1.03}]},
		{GrayLevel[1.],Hue[.95],Text["d",{1.1,-1.1}]}};
		
EarthTime[t_Real] :=				{GrayLevel[1.],Hue[.7],Text["Earth Time: "<>
				ToString[NumberForm[t,{5,2}]]<>" y",{-.95,-1.1},{-1,0}]};
ShipTime[t_Real]	:=
				{GrayLevel[1.],Hue[.7],Text["Ship Time: "<>
					ToString[NumberForm[t,{5,2}]]<>" y",{.05,-1.1},{-1,0}]};
VMeter[v_Real]		:=				{GrayLevel[1.],Hue[.95],Rectangle[{-1.14,-1.0},{-1.06,2*v-1.}]};
DMeter[d_Real]		:=
				{GrayLevel[1.],Hue[.95],Rectangle[{1.06,-1.0},{1.14,2*d/dTrip-1.}]};
Title[sname_String,dist_Real] :=	
				{GrayLevel[1.],Hue[.7],
				Text["Destination: "<>sname<>" "<>ToString[NumberForm[dist,{4,1}]]<>
							" ly",{0.,1.1}]};

Bring in star catalogue

While[
(StarFile = InputString["Enter name of star catalogue:\n"<>
				"bsc4.m  (stars up to 4th mag.)\n"<>
				"bsc5.m  (stars up to 5th mag.)\n"<>
				"bsc6.m  (stars up to 6th mag.)\n"];
strm = OpenRead[StarFile])===$Failed,];
StarCatalogue = ReadList[strm,{Word,Word,Number,Number,Number,Number},
WordSeparators-> {"\t"},RecordSeparators->{"\n"}];
Close[strm];

Star[0] = {"Sun","",0.0,90.,-26.74,begpt};
Star[nn_ ]/; nn >= 0 := StarCatalogue[[nn]] ;
phi[nn_ ]/;nn >= 0 := N[Star[nn][[3]] Degree] ;
tht[nn_ ]/;nn >= 0 := N[Star[nn][[4]] Degree] ;
nStars = Length[StarCatalogue];
Do[uzStar[k] =  Sin[tht[k]]{Cos[phi[k]],Sin[phi[k]],0} +
				Cos[tht[k]]{0,0,1},
	{k,0,nStars}];
Print["Catalogue ",StarFile," read with ",nStars," stars."];

(*  Defining our trip  *)

The lorentz boost at constant acceleration g*g0 from M.T.W. p166 is:

Clear[DefineTrip];
DefineTrip[n_Integer,g_Real,gmx_Real,plot_:True] := 		Block[{tp = t0/g,tpp, dcrit = 2*tp*(gmx-1)},
			{sname,constel,ra,codec,smag,invdly} = Star[n];
			dly = 1/invdly;
			Print["Setting up trip to ",sname," in ",constel];
			Print["Distance: ",NumberForm[dly,{5,2}]," light-years."];
			BackgroundFrame = Append[MovieFrame,Title[sname,dly]];
			dTrip = dly - endpt;
			If[dcrit< dTrip,
				(
				tauAcc		= tp ArcCosh[gmx];
				tauCoast	= (dly - 2*tp*(gmx-1))/Sqrt[gmx^2-1];
				),
				(
				tauAcc 	=  tp ArcCosh[(dTrip/(2*tp) + 1)];
				tauCoast	= 0.;
				)
				];
			gammaMax = Cosh[tauAcc/tp];
			betaMax	= Tanh[tauAcc/tp];
			tauTrip	= 2*tauAcc + tauCoast;
			Print["Ship time for trip: ",NumberForm[tauTrip,{5,2}]];
			timeAcc		= tp*Sinh[tauAcc/tp];
			timeTrip		= 2*timeAcc + gammaMax*tauCoast;
			Print["Earth time for trip: ",NumberForm[timeTrip, {5,2}]];
			rangeAcc		= tp*(Cosh[tauAcc/tp] -1);
			thetaTrip  =  codec Degree // N;
			phiTrip     =  ra Degree // N;
			uzTrip =					Sin[thetaTrip]{Cos[phiTrip],Sin[phiTrip],0} +
					Cos[thetaTrip]{0,0,1};
			uxTrip = {-Sin[phiTrip],Cos[phiTrip],0};
			uyTrip = Cos[thetaTrip]{Cos[phiTrip],Sin[phiTrip],0} -
					Sin[thetaTrip]{0,0,1};
			uTrip = {uxTrip,uyTrip,uzTrip};
			uTripInv = Transpose[uTrip];
(* Trip functions  *)
			Clear[range,time,gamma];
			range[tau_]=						Which[
							tau < 0,
									0.,							tau <= tauAcc,
									tpp*(Cosh[tau/tpp]-1),
							tau <= tauAcc + tauCoast,									rangeAcc + Sqrt[gammaMax^2-1](tau-tauAcc),
	 						tau <= tauTrip,
	 								dTrip - tpp*(Cosh[(tauTrip-tau)/tpp]-1),
	 						True,
	 								dTrip
	 						] /. tpp -> tp;
			time[tau_]= 					Which[
							tau < 0,
									tau ,
							tau <= tauAcc,
									tpp*Sinh[tau/tpp],
							tau <= tauAcc + tauCoast,
									timeAcc + gammaMax*(tau-TauAcc),
							tau <= tauTrip,
									timeTrip - tpp*Sinh[(tauTrip-tau)/tpp],
							True,
									timeTrip +(tau - tauTrip)
							] /. tpp -> tp;
			gamma[tau_]= 					Which[
							tau < 0,0.,
							tau <= tauAcc,Cosh[tau/tpp] ,
							tau <= tauAcc + tauCoast, gammaMax,
							tau <= tauTrip,Cosh[(tauTrip-tau)/tpp],
							True,0.]/. tpp -> tp;
							
			beta[tau_]:= Sqrt[1-1/gamma[tau]^2];
			If[plot,
			Plot[time[t],{t,0.,tauTrip},Frame->True,
				FrameLabel->{"Proper Time (yrs)","Observer Time (yrs)"} ];
			Plot[range[t],{t,0.,tauTrip},Frame->True,
				FrameLabel->{"Proper Time (yrs)","Distance (light-yrs)"} ];
			Plot[beta[t],{t,0.,tauTrip},Frame->True,PlotRange->{{0.,tauTrip},{0.,1.}},
				FrameLabel->{"Proper Time (yrs)","Beta"} ];
			Plot[gamma[t],{t,0.,tauTrip},Frame->True,
					PlotRange->{{0.,tauTrip},{0.,gammaMax}},
					FrameLabel->{"Proper Time (yrs)","Gamma"} ]
				];
		];

(* Vector functions *)
Mag[v_List] := Sqrt[v.v];

(*  Camera direction set up *)
thCamOld = 0.;
phCamOld = 0.;
uxCam = {1.,0.,0.};
uyCam = {0.,1.,0.};
uzCam = {0.,0.,1.};

Clear[SetCamera];
SetCamera[n_Integer,fnum_:1.] :=			Block[{sname,constel,magstar,invdstar},
			{sname,constel,phcam,thcam,magstar,invdstar} = Star[n];
				phcamr = phcam Degree // N;
				thcamr = thcam Degree // N;
				fCam = fnum;
				If[thcam =!=thCamOld || phcam =!= phCamOld,
					{uxCam,uyCam,uzCam} =					{{Sin[phcamr],-Cos[phcamr],0.},
					{-Cos[thcamr]Cos[phcamr],-Cos[thcamr]Sin[phcamr],Sin[thcamr]},
					{Sin[thcamr]Cos[phcamr],Sin[thcamr]Sin[phcamr],Cos[thcamr]}};
					thCamOld = thcam;
					phCamOld = phcam]
				];

(* Star plotting  *)
Color[mag_,dop_] :=				Block[{br,sat},
					br = Min[1., .7/(1. + 10.^((mag-2)/2.)) +.3];
					ld = Log[5.,dop];
					col = If[ld>0.,
							Hue[.60,Min[1.,ld],br],
							Hue[1.0,Min[1.,-ld],br]
							];
		Return[col]		
				];
		
Size[mag_] := AbsolutePointSize[2];

Clear[StarImage];
StarImage[mag_,dop_,{x_,y_}] :=									{Size[mag],Color[mag,dop],Point[{x,y}]};
			
Clear[PlotStar];
PlotStar[n_Integer] :=
		Block[ {},				(* {sname,constel,phstar,thstar,magstar,dstar,
				bri,obs,kstar,magobs,vecstar,
				gam,bb,Kstar,Kstarp,kk,doppler} *)
					{sname,constel,phstar,thstar,magstar,invdstar} = Star[n];
			 		
				  kstar	=	If[invdstar == 0,
				  						(magobs = magstar;
				  							uzStar[n]),
				  						(vecstar = uzStar[n]/invdstar - rangeNow*uzTrip;
				  							magobs = magstar + 5.*Log[10.,(invdstar*Mag[vecstar])];
				  							kstar = vecstar/Mag[vecstar]
				  							)
				  						];
				  gam = gammaNow;
				  bb	= betaNow;
				  Kstar = Append[uTrip.kstar,1.]//Chop;
				  Kstarp =
				  		{{1,0,0,0},{0,1,0,0},{0,0,gam,bb*gam},{0,0,bb*gam,gam}}.
				  			Kstar//Chop;
				  doppler = Kstarp[[4]];
				  	kkp = uTripInv.(Drop[Kstarp,-1]/doppler)//Chop;
				  	kkx =  kkp.uxCam // Chop;
				  	kky =  kkp.uyCam // Chop;
				  	kkz =  kkp.uzCam // Chop;
				  	If[kkz <= 0., Return[$Failed]];
				  	xcam = fCam*kkx/kkz;
				  	ycam = fCam*kky/kkz;
				  	If[Abs[xcam] > 1. || Abs[ycam] > 1.,Return[$Failed]];
				  	rslt = StarImage[magobs,doppler,{xcam,ycam}];
				  	Return[rslt]
				  		]

(*  Frame definition  *)
StarFrame[tau_Real]:=				Block[{tt=tau,frame,stars,nn},
						rangeNow = range[tt];
						gammaNow = gamma[tt];
						betaNow		= Sqrt[Max[gammaNow^2-1,0.]]/gammaNow;
						timeNow		= time[tt];
						frame = 								Join[BackgroundFrame,
									{EarthTime[timeNow],ShipTime[tt],
									VMeter[betaNow],DMeter[rangeNow]}];
						stars = {};
						Do[
							(
							result=PlotStar[nn];
							If[result=!=$Failed,AppendTo[stars,result]]
							),
						{nn,0,nStars}];
						frame = Graphics[Append[frame,{GrayLevel[1.],Reverse[stars]}]];
						Return[frame]
						];

SingleFrame[tau_Real] :=				Show[StarFrame[tau],
					AspectRatio->Automatic,
					PlotRange->{{-1.2,1.2},{-1.2,1.2}}];
													
MakeMovie[nframes_Integer] :=
	Do[SingleFrame[ttt],{ttt,0.,tauTrip,tauTrip/nframes }];

EndPackage[];