{~ Part of the Map_BSC programme, split from it to enable debugging ~} unit {map_bs_s} map_bsts {replacement of map_bs_s, running in text mode} {it should not be edited, as it is always generated by a change file from the file it replaces} ; {used just by the Map_BSC programme} interface uses dos; const rem_char=['&','%','#']; step_plan:double=1; par_file_Name:Pathstr=''; Morning_or_evening:char=' '; {M for Morning, E for Evening, B for Both} Faint_abs_cmag:integer=-2000; Sun_depth:double=12; Write_time:boolean=true; {Show time by each position} Stereo_Dist_ftEarth:double=0; nst: string=''; NoS: longint=0; WtTimeMod:byte=1; WtTimeModMoon:byte=1; Faintness_lim_SolSyst:word=800; OutGas:boolean=false; MoreTimesHappened:boolean=false; ShowDistances:boolean=false; ShowGeocentricMoonToo:boolean=true; ShowProgress:boolean=true; PlanFullMode:boolean=false; {default empty planets, parameter PLF overrides it} var Outside:word; Tyc_str:string[25]; {just for Progress report for Tycho} remark,remark_line_file,LastReport:boolean; par_line,command_line,name:string; disp_times:boolean; eph_time,OSM,OJD,FirstJD,FirstSM,LastJD,LastSM, bg,rg,Fi_true,Sol_rad,step: real; MoE:char; par_file:text; xnum,report_h,V2RV:word; SigStep:shortint; procedure posi_elm; procedure plan_moon(Jd_p,sm_p:double;n_plp:byte;Moon_only:boolean); Procedure Display_time(FirstDesc:boolean;WriteTime:boolean; JD_mom:double;shift_right:integer); procedure progress_report; {informs on the plotted/read number of stars} function ShowAsStar(Az:integer;Empt:boolean;rg:double;lab_st:string):boolean; Procedure NewTime; procedure FirstLast; implementation uses crt,{graph} graph_tx, {J.H.'s text replacement of graph} str_num,{params_g} params {replacement of params_g, running in text mode} ,{graph_m8} graph_mt {replacement of graph_m8, to avoid screen graphics} {it should not be edited, as it is always generated by a change file from the file it replaces} ,angles_o,un_str_m, solar_ut,planets,{map_bs_p} map_bstp {replacement of map_bs_p, running in text mode} {it should not be edited, as it is always generated by a change file from the file it replaces} ,precesse,elm; {TYPE V3 = array[1..3] of double;} const DaysOfWeek:boolean=true; First_plan:array[0..10] of boolean= (true,true,true,true,true,true,true,true,true,true,true); Shift_Time:integer=83; var Day:string; YYYYold,MMold:word; NDold:integer; procedure progress_report; {informs on the plotted/read number of stars} begin if Monoscopic and ShowProgress then begin setcolor(0); moveto(xnum,report_h); outtext(nst); setcolor(wrtcolour); moveto(xnum,report_h); nst:=SI(4,NoPlotted)+' / '+SI(6,NoS); if LastReport then begin setcolor(0); moveto(5,report_h); outtext('plotted/read, free mem.: '); setcolor(wrtcolour); moveto(5,report_h); outtext('plotted/read: '); xnum:=getx; end else nst:=nst+SI(8,MaxAvail); if Tycho and (not LastReport) then nst:=nst+'; '+Tyc_str+SI(5,Outside)+' (C for continue)'; outtext(nst); end; end; Procedure NewTime; begin if Morning_or_evening=' ' then begin eph_time:=eph_time+step*SigStep; JD:=JD+step*SigStep; SM:=Bas_Int(SM+24.0657098230*step*SigStep,24); end else begin if morning_or_evening='B' then if MoE='M' then MoE:='E' else begin eph_time:=eph_time+step*SigStep; MoE:='M' end else eph_time:=eph_time+step*SigStep; if not At_height(MoE,-Sun_depth,eph_time) then halt; JD:=eph_time+JD_MJD; Name:=DoW(JD,t_UT,_English); DD:=ND; Get_JD; if surface =cylinder then ST:=SM; end end; procedure FirstLast; begin if JD>LastJD then begin LastJD:=JD; LastSM:=SM; MoreTimesHappened:=true; end else if JD spoil vertical centering} else Day:='('+copy(Day,1,3)+')'; Day:=__Dt2StStr(YYYY,MM,ND,_ANSI_Dt_Str)+Day; if FirstDesc then begin YYYYold:=YYYY; MMold:=MM; NDold:=ND end else if YYYYold=YYYY then begin Day:=copy(Day,5,length(Day)-4); if MMold=MM then begin Day:=copy(Day,4,length(Day)-3); if NDold=ND then Day:=' ' else NDold:=ND; end else begin MMold:=MM; NDold:=ND; end end else begin YYYYold:=YYYY; MMold:=MM; NDold:=ND; end; end else begin Day:=''; Time:=''; name:=__cvtstr(copy(name,2,length(name)-1),_rem_white_str); end; gm_outxtxy(x,y,'',Day+Time+name,'','l',round_rot(A),shift_right+16,0) end; function ShowAsStar(Az:integer;Empt:boolean;rg:double;lab_st:string):boolean; var VTV:integer; begin if V>Faintness_lim_solsyst then ShowAsStar:=false else begin VTV:=VT(V); V2RV:=VdmagF2R(VTV); ShowAsStar:=true; if stereoscopic then begin Parallax:=round(stereo_coef_plan/(rg-Stereo_Dist_ftEarth)); x:=(x div 16)*16 + 11; {for getting smooth one-pixel screen parallax step, but coarser stellar positions in PostScript file as well} x_left:=x+Parallax; if StereoMirror then x_left:= StereoHalfX - x_left; gm_star(x_left, y,VTV,V2RV,Az,Empt,''); if lab_st<>'' then gm_outxy(x_left,y,lab_st,'c',round_rot(A)); if StereoColor then begin if Gray16_stars then begin StereoRight:=true; setcolor(wrtcolour shl TwoOrFour) end else SetColor(right_eye {lightgreen/blue} ); if GM_mode=3 then gm_verbatim('/pc {pcr} def pc '); gm_star_full_mode:=OrPut; {to prevent biting into the left star; red transparent for black would be better, of course, then even faint the green/blue stars would be visible within the bright ones} end; x_left:=x+stereo_shift-Parallax; gm_star(x_left,y,VTV,V2RV,Az,Empt,''); if lab_st<>'' then gm_outxy(x_left,y,lab_st,'c',round_rot(A)); if StereoColor then begin if Gray16_stars then begin StereoRight:=false; setcolor(wrtcolour) end else SetColor(LightRed); if GM_mode=3 then gm_verbatim('/pc {pcl} def pc '); gm_star_full_mode:=star_full_mode_ori; end; end else begin gm_star(x,y,VTV,V2RV,Az,Empt,''); if lab_st<>'' then if not PlanFullMode then gm_outxy(x,y,lab_st,'c',round_rot(A)) else gm_outxtxy(x,y,'',lab_st,'','r',round_rot(A), -V2RV-16,0); end; end; end; procedure posi_elm; const hp10= cl+ '$@'+cl+ 'should contain (and parameters given elsewhere, preceded by $, can change):'+cl+ cl+ '{ {A|Q}# AN# I# E# P{D|M#} P{N|L}#'+cl+ ' | ![:] }'+cl2; hp15= ' ![:] will read elements from the line containing '+cl+ ' using text(!) file (default xephem.dat, mpc format is,'+cl+ ' an alternative, in this case >mpc< should be contained in the file name'+cl2; hp20= ' A is the large semiaxis'+cl+ ' AN longitude of ascending node'+cl+ ' I inclination'+cl+ ' E numerical excentricity (default 0)'+cl+ ' PD perihelion date'+cl+ ' PM perihelion Modified Julian Date'+cl; hp25= ' (moment of perihelion passage can be replaced by:'+cl+ ' MA mean anomaly, AD that mean Anomaly Date and ev. DM daily mean motion)'+cl; hp30= ' PN elongation of perihelion from ascending node'+cl+ ' PL perihelion longitude (= PN + AN)'+cl+ ' Q perihelion distance from the Sun'+cl2; hp40= 'It may contain also moments and way the object should be shown:'+cl+ ' dDD[.MM[.YY]] tHH[.MM[.SS]] UT# FA#[:#] CE[#] B# N# NT S{#|#h|#m} M'+cl+ ' CD# ST#'+cl; hp41= ' t{A#:#:... | B# | D#[:#[...]] | I[#] | M# |'+cl+ ' O[#] | P#:#:# | T#[:#]~#:#~...| W#[:#[:#]] }'+cl+ ' { ? | h }' +cl2; hp45= ' D and T are initially taken from the main programme'+cl+ ' UT is the t-UTC difference / 1 h (November-March 1, then 2)'+cl; hp50= ' FA Absolute Faintness / 1 cmag (those at 1 AU from both Earth and Sun,'+cl+ ' at zero "phase angle" (in fact, Sun-object-Earth <> 0 at this position)'+cl+ ' [: (def. 400)]'+cl; hp55= ' B positions before the given moment (default 0)'+cl; hp60= ' N next positions to be computed (default 0)'+cl+ ' NT[-] do not (-: do) show times (and dates) by each position'+cl+ ' S step to further/preceding time (default 1 day, S#h for time in hours)'+cl+ ' M string to be shown at each place (default +)'+cl; hp65= ' W module for writing time (default as for planets)'+cl+ ' [: (default 10* the first one)'+cl+ ' [: (default the preceding one /5)]]'+cl2; hp70= 'The remaining parameters concern comets only:'+cl; hp80= ' CE cometary exponent (amount of diffused sunlight is proportional to'+cl+ ' distance from Sun **CE, default 4, solid bodies have 2); phase angle'+cl+ ' dependence of faintness is neglected when the parameter is given'+cl; hp90= ' CD Time of negligible solar acceleration of particles when blown by gas'+cl+ ' in coma (default 2 days)'+cl+ ' ST Split Time of large dust particles (default no splitting)'+cl2; hp100= 'parameters beginning with T (Tail; finish its computation by pressing >S<):'+cl; hp110= ' tA set of solar repulsion accelerations of cometary particles, as ratios'+cl+ ' to gravitational solar attraction, maximum 6 values separated by :'+cl+ ' (default 0.5:1:2:3:4:200), less then 6 values imply less particles;'+cl; hp111= ' repulsion may be appended by a value differing from that given by CD'+cl+ ' (delimiter is ", i.e., last item may read 200"4 instead of just 200)'+cl; hp120= ' tB number of parent body positions before the given moment for which'+cl+ ' cometary particles should be generated (default = B)'+cl+ ' tD ejection velocity of particles at 1 AU / (m/s) (default 400);'+cl+ ' more values (separated by :) apply to last tA values'+cl; hp130= ' tI interrupt computation each # step (1 when no # is given, pressing >R<'+cl+ ' begins continuous run of programme without interruptions)'+cl+ ' tM show just tails mod # (default=tB, just the tail in the given moment'+cl+ ' is shown)'+cl; hp140= ' tO show tail Only (do not mark position of the parent body); if time is' +cl+ ' to be written (see $W parameter), shift it # millimeters right'+cl; hp145= ' tP (Prograde) pole of the rotation of the parent body and longitude of axis'+cl+ ' of ejection measured from the subsolar meridian (RA:decl:angle)/1 degree'+cl; hp150= ' tT Tilts of ejection (from its axis) / 1 degree'+cl+ ' and number of sections of the cone given as tilt[:sect.]~tilt[:sect.]~...'+cl+ ' (default tilt 0, sections 6 for tilt>0)'+cl; hp160= ' tW velocity of the solar wind / 1 km.s (default 400)'+cl; hp170= ' ? will display this help (i.e., $? from the command line)'+cl+ ' h will append it to the Map_BSC.hlp file'+cl2; hp180= 'Each $@ represents one body.'+cl+ ' The same holds for $!: -- an alternative for a $@'+cl+ ' as far as the necessary ephemeris parameters are concerned.'+cl; hp190= 'If more than one is given in a @'+cl+ ' or on the command line,'+cl+ ' only those parameters are valid for the current ${@|!}file,'+cl+ ' which are given before the following ${@|!}file.'+cl; Next:word=0; Before:word=0; {Kind:byte=3;} {O_Name:string='object';} Symb:string[10]='+'; NewMJD:boolean=false; {serves for possibility to plot [comets/minor planets] of different years (or of another times in general) into the same map} Body_No:byte=0; Elm_file_No:byte=0; TailMod:byte=0; TailBefore:byte=0; n_Acc:byte=6; n_tilts:byte=1; {a_beh:double=0;} PartAcc:array[1..6] of double = (0.5,1,2,3,4,200); TailWait:byte=0; coma_days:array[1..6] of double= (2,2,2,2,2,2); Rot_angle:double=0; vDust:array[1..6] of double=(400,400,400,400,400,400); {m/s} vSolWind:double=400; {km/s} SplitTimeTrue:double=0; {days} PoleGiven:boolean=false; old_elm_file_Name:Pathstr=''; TailOnly:boolean=false; type TiltRec= record tiltr:double; n_dir:byte; end; var PartTilt:array [1..4] of TiltRec; const Tilts:boolean=false; var xe,ye,ze,l,b,drp:real; mp,dp:word; yp:integer; i:byte; elm_file_name:Pathstr; TailPole:V3; Found_body:boolean; Write_time_body,Write_time_tail:boolean; WtTimeMod_body,WtTimeMod_tail,WtTimeMod_tail_g:byte; Faint_add:word; procedure Help_p(s:string); var ch:char; begin if gm_open then gm_close; if s<>'' then writeln(s); write(hp10,hp15,h0); ch:=readkey; if UpCase(ch)= 'Q' then halt; write(h0clear,hp20,hp25,hp30,h0); ch:=readkey; if UpCase(ch)= 'Q' then halt; write(h0clear,hp40,hp41,hp45,hp50,hp55,hp60,hp65,hp70,h0); ch:=readkey; if UpCase(ch)= 'Q' then halt; write(h0clear,hp80,hp90,hp100,hp110,hp111,h0); ch:=readkey; if UpCase(ch)= 'Q' then halt; write(h0clear,hp120,hp130,hp140,hp145,hp150,hp160,h0last); ch:=readkey; if UpCase(ch)= 'Q' then halt; write(h0clear,hp170,hp180,hp190); halt end; procedure help_p_file; var h_file:text; begin if rt_append('',h_file,'Map_BSC.hlp') then begin write(h_file, hp10,hp15,hp20,hp25,hp30, hp40,hp41,hp45,hp50,hp55,hp60,hp65,hp70, hp80,hp90,hp100,hp110,hp111, hp120,hp130,hp140,hp145,hp150,hp160, hp170,hp180,hp190); close(h_file); gm_close; halt end else help_p(''); end; procedure one_p_param(inside_elm_file:boolean); var j_ic,j_aux:byte; aux_st:string[10]; ParOrgCase:string; begin if (UpCase(Par[1])='M') and (length(Par)>1) then Symb:=copy(Par,2,length(Par)-1) else begin ParOrgCase:=Par; Par:=__cvtstr(Par,_To_UpCase_Str); case Par[1] of '?' : Help_p(''); 'H' : Help_p_file; '!': if inside_elm_file then begin if ItemCountD(':',Par)>1 then ElmCatNameDefault:=ItemStrD(':',2,ParOrgCase); if pos('mpc',ElmCatNameDefault)>0 then ElmCatNameKind:=2; if not get_orbital_elements( ElmCatNameKind, ElmCatNameDefault, ItemStrD(':',1,copy(ParOrgCase,2,length(par)-1))) then help_p(ItemStrD(':',1,copy(ParOrgCase,2,length(par)-1))+' not found'); Faint_abs_cmag:=roundup(Faint_abs*100); end; 'A': case Par[2] of 'N': asc_node:=deg_to_rad(ss2r(Par,3)); 'D': begin {YY:MM:DD.dddd decoding} __St2RDtStr(copy(Par,3,length(Par)-1),YP,MP,DP,DRP,_ANSI_Dt_Str); anomaly_date:=JD0(YP,MP,trunc(DRP))+frac(DRP)-JD_MJD end; else semiaxis_large:=Par2r(par); end; 'B': if length(par)>1 then Before:=ss2i(par,2); 'C': if length(par)>1 then case Par[2] of 'E': begin if length(Par)>2 then Comet_exponent:=ss2r(Par,3) else Comet_exponent:=4; Comet:=true; end; 'D': if length(Par)>2 then begin Coma_days[1]:=ss2r(Par,3); for j_ic:=2 to 5 do Coma_days[j_ic]:=Coma_Days[1]; end; end; 'D': if length(Par)>1 then if Par[2]='M' then begin if length(Par)>2 then daily_motion:=ss2r(Par,3); end else if not Par_Date then help_p('') else NewMJD:=true; 'E': if length(par)>1 then excentricity:=Par2r(par); 'F': if length(par)>2 then if Par[2]='A' then begin Symb:=''; Faint_abs_cmag:=ss2i(ItemStrD(':',1,Par),3); if ItemCountD(':',Par)>1 then Faint_add:=I_S(ItemStrD(':',2,Par)) else Faint_Add:=400; end; 'M': if length(Par)>2 then if Par[2]='A' then mean_anomaly:=ss2r(Par,3); 'N': if length(par)>1 then if Par[2]='T' then if length(par)>2 then if par[3] in ['-','F'] then begin if WtTimeMod_body>0 then Write_time_body:=true end else Write_time_body:=false else Write_time_body:=false else Next:=Par2i(Par); 'I': if length(par)>1 then inclin:=deg_to_rad(Par2r(par)); 'P': if length(par)>1 then case Par[2] of 'N': perih_from_asc_node:= deg_to_rad(ss2r(Par,3)); 'L': perih_comp_length:= deg_to_rad(ss2r(Par,3)); 'D': begin {YY:MM:DD.dddd decoding} __St2RDtStr(copy(Par,3,length(Par)-1),YP,MP,DP,DRP,_ANSI_Dt_Str); perih_pass:=JD0(YP,MP,trunc(DRP))+frac(DRP)-JD_MJD end; 'M': perih_pass:=ss2r(Par,3); end; 'S':if length(par)>1 then case par[2] of 'T': SplitTimeTrue:=ss2r(par,3); else begin if pos('H',Par)>0 then Step:=par2r(copy(Par,1,pos('H',Par)-1))/24 else if pos('M',Par)>0 then Step:=par2r(copy(Par,1,pos('M',Par)-1))/24/60 else Step:=par2r(Par); end; end; 'Q': if length(par)>1 then perihel_dist:=par2r(Par); 'T': if length(par)>1 then case par[2] of 'A': if length(par)>2 then begin par:=copy(par,3,length(par)-2); n_Acc:=ItemCountD(':',par); if n_Acc>6 then n_Acc:=6; for j_ic:=1 to n_Acc do begin aux_st:=ItemStrD(':',j_ic,par); if ItemCountD('"',aux_st)>0 then begin Coma_days[j_ic]:=R_S(ItemStrD('"',2,aux_st)); aux_st:=ItemStrD('"',1,aux_st); end; PartAcc[j_ic]:=R_S(aux_st); end; end; 'B': TailBefore:=ss2i(par,3); 'D': if length(par)>2 then begin par:=copy(par,3,length(par)-2); j_ic:=ItemCountD(':',par); if j_ic>n_Acc then j_ic:=n_Acc; for j_aux:= j_ic downto 1 do vDust[n_Acc+j_aux-j_ic]:=R_S(ItemStrD(':',j_aux,par)); if j_ic2 then TailWait:=ss2i(par,3) else TailWait:=1; 'M': TailMod:=ss2i(par,3); 'O': begin TailOnly:=true; if (TailBefore=0) and (Before>0) then TailBefore:=Before; Before:=0; if length(par)>2 then Shift_Time:=round(ss2r(par,3)/gm_un2mm); end; 'P': if length(par)>2 then begin par:=copy(par,3,length(par)-2); if ItemCountD(':',par) = 3 then begin l:=deg_to_rad(R_S(ItemStrD(':',1,par))); b:=deg_to_rad(R_S(ItemStrD(':',2,par))); ra_d2l_b(epsi_rad,l,b); P3s_vector(l,b,TailPole); Rot_angle:=deg_to_rad(R_S(ItemStrD(':',3,par))); PoleGiven:=true; end; end; 'T': if length(par)>2 then begin par:=copy(par,3,length(par)-2); n_Tilts:=ItemCountD('~',par); if n_Tilts>4 then n_tilts:=4; for j_ic:=1 to n_Tilts do with PartTilt[j_ic] do begin aux_st:=ItemStrD('~',j_ic,par); tiltr:=R_S(ItemStrD(':',1,aux_st)); if tiltr>0 then begin tiltr:=deg_to_rad(tiltr); Tilts:=true; end; if ItemCountD(':',aux_st)>1 then n_dir:=I_S(ItemStrD(':',2,aux_st)) else if tiltr=0 then n_dir:=1 else n_dir:=6; if n_dir=0 then inc(n_dir); end; end; 'W': if length(par)>2 then vSolWind:=ss2r(par,3); else if not Par_Time then help_p('') else NewMJD:=true; end; 'U': if not Par_UT then help_p(''); 'W': if length(par)>1 then begin WtTimeMod_body:=par2i(ItemStrD(':',1,Par)); if ItemCountD(':',Par)>1 then begin WtTimeMod_tail:=I_S(ItemStrD(':',2,Par)); if ItemCountD(':',Par)>2 then WtTimeMod_tail_g:=I_S(ItemStrD(':',3,Par)) else WtTimeMod_tail_g:=WtTimeMod_tail div 5; end else if WtTimeMod_body<11 then begin WtTimeMod_tail:=WtTimeMod_body * 10; WtTimeMod_tail_g:=WtTimeMod_body * 2; end else begin WtTimeMod_tail:=0; WtTimeMod_tail_g:=0; end; Write_time_body:= WtTimeMod_body >0; Write_time_tail:= (WtTimeMod_tail >0) and (WtTimeMod_tail_g>0); end; end; {of case} end; end; procedure Init_next_bef; begin PartTilt[1].tiltr:=0; PartTilt[1].n_dir:=1; step:=step_plan; Write_time_body:=Write_time; Write_time_tail:=Write_time_body; WtTimeMod_body:=WtTimeMod; WtTimeMod_tail:=WtTimeMod_body * 10; WtTimeMod_tail_g:=WtTimeMod_body * 2; end; procedure Read_elm_file; var jjj,ic:byte; elm_remark:boolean; elm_file:text; elm_line:string; begin Init_next_bef; if not rt_reset_suf(Loc_Obj,elm_file,elm_file_name,'elm') then help_p('@ file >>'+elm_file_name+'<< not found'); while not eof(elm_file) do {processing the $@} begin readln(elm_file,elm_line); elm_remark:=false; ic:=itemcount(elm_line); for jjj:= 1 to ic do begin Par:=itemstr(jjj,elm_line); if Par[1] in rem_char then elm_remark:=true; if not elm_remark then one_p_param(true); end; end; close(elm_file); elm_remark:=false; end; procedure do_body; const BeforeComputing:boolean=false; p_start_po:pointer=nil; var N_other:word; MJD_start:double; DoTail:boolean; debug_ch:char; heading:string[50]; ShowPosition{,StopTail}:boolean; procedure NextBefore; type Particle=record xp: V3; vp: V3; old: byte; {Number of steps (default = 1 d) from the given ephemeris time} kin: byte; central:boolean; fa: integer; nextP: pointer; end; Position=record eph_time_po:double; ST_po:double; next_po:pointer; end; cy3=0..2; const j_ori:cy3=0; var j_old:byte; jjj,j_tilt,j_dir,SplitTime:byte; j:integer; jcoord:c3d; jori_m1,jori_m2:cy3; VAcc,VVel,VTrajPole,VVelEj:V3; GSol,vSolWind2,r3,v_ej,eph_time_comp:double; fa_ini:integer; p_part,p_part_new,p_part_ins:^particle; p_posi,p_posi_new:^position; p_start:pointer; {to the first old particle} xo:array[0..2] of V3; ro:array[0..2] of double; MVel,MRot:M3; procedure ShowTailParticle; begin with p_part^ do begin c3_to_Pd(xp[1],xp[2],xp[3],r,l,b); Hel_to_Geo(eph_time_comp+JD_MJD+tt_ut,l,b,r,lg,bg,rg); LB_RD(lg,bg,ra_d,de_d); RADe2xy; if aux_boo then begin V:=fa+RoundUp(500*ln(r*rg)/l10); if Elevate_by_refraction then inc(V,round(ZenExt*Air_Mass)); if ShowAsStar(0,false,rg,'') then if Write_time_tail and central and (old>0) and (((old mod (WtTimeMod_tail) = 0) and (old > WtTimeMod_tail) and (kin in [1..n_Acc-1])) or ((kin in [n_Acc,0]) and (old mod WtTimeMod_tail_g = 0) and (old > WtTimeMod_tail_g))) then begin if ShowDistances then name:='d;'+SRx(1,r)+'AU;'+SRx(1,rg)+'AU' else name:=''; gm_outxtxy(x,y,'', __cvtstr(SI(3,RoundUp(old*step)),_rem_white_str)+name,'', 'c',round_rot(A),0,0); end; end; end; end; {of procedure ShowTailParticle, within NextBefore} procedure ShiftOldParticle; var jcoord:c3d; VAccL:v3; procedure add_splitted_particle(dust:boolean); var jc:c3d; begin with p_part^ do begin new(p_part_ins); p_part_ins^.nextP:=nextP; if dust then begin p_part_ins^.kin:=kin; p_part_ins^.old:=old; end else begin p_part_ins^.kin:=n_Acc; p_part_ins^.old:=0; end; p_part_ins^.fa:=fa; p_part_ins^.central:=central; if p_part_ins^.kin vSolWind2 then kin:=0; if kin>0 then for jcoord:= 1 to 3 do begin if step*old>coma_days[kin] then {no repulsive accelaration yet} VAccL[jcoord]:=-xp[jcoord]*(1-PartAcc[kin])*GSol/r3 else {repulsive acceleration already} VAccL[jcoord]:=-xp[jcoord]* GSol/r3; xp[jcoord]:=xp[jcoord]+vp[jcoord]+VAccL[jcoord]/2; vp[jcoord]:=vp[jcoord]+VAccL[jcoord] end else for jcoord:= 1 to 3 do xp[jcoord]:=xp[jcoord]+vp[jcoord]; inc(old); if (j-j_old) mod TailMod = 0 then ShowTailParticle; if (SplitTime>0) and (kin<3) and (old mod SplitTime = 0) and (kin0) then begin MJD:=MJD_start; eph_time:=MJD; if (surface=cylinder ) and (Morning_or_evening<>' ') then ST:=OSM; end; if DoTail then begin new(p_part); p_part^.nextP:=nil; p_start:=p_part; GSol:=3E-4*sqr(step); vSolWind2:=sqr(vSolWind*86400000*step/AU ); for jjj:=1 to n_Acc do vDust[jjj]:=vDust[jjj]*86400*step/(AU); SplitTime:=round(SplitTimeTrue/step); P3s_Vector(asc_node-pi/2,pi/2-inclin,VTrajPole); end; if BeforeComputing then SigStep:=-1 else SigStep:= 1; for j:=0 to N_other do begin object_hel( semiaxis_large, perihel_dist, excentricity, inclin, asc_node, perih_from_asc_node, perih_pass, eph_time+tt_ut, xe,ye,ze,r); ShowPosition:= ((j>0) or (not BeforeComputing) or (Next=0)) {base date just once} and ((not BeforeComputing) or (j<=Before)) {not just tail production} and (not TailOnly); {not skipping circle of position} if ShowPosition or ((Write_time_body or ShowDistances) and (j mod WtTimeMod_body = 0)) then begin c3_to_Pd(xe,ye,ze,r,l,b); Hel_to_Geo(eph_time+JD_MJD+tt_ut,l,b,r,lg,bg,rg); if ShowDistances then begin name:=SRx(1,r)+'AU;'+SRx(1,rg)+'AU'; if Write_time_body then name:=';'+name; end else name:=''; LB_RD(lg,bg,ra_d,de_d); RADe2xy; if aux_boo then begin if ShowPosition then begin if Symb<>'' then begin gm_outxy(x,y,Symb,'c',round_rot(A)); V2RV:=textwidth(symb)*16; end else begin if Comet then V:=Faint_abs_cmag+RoundUp(500*((Comet_exponent/2)*ln(r)+ln(rg))/l10) else V:=Faint_abs_cmag+round(100* planet_faintness(0,object_phase_angleD(l,b,lg,bg),r,rg)); if Elevate_by_refraction then inc(V,round(ZenExt*Air_Mass)); ShowPosition:=ShowAsStar(round_rot(A),true,rg,''); end; if ShowPosition and (Write_time_body or ShowDistances) and (j mod WtTimeMod_body = 0) then Display_time(j=0,Write_time_body, eph_time+JD_MJD,V2RV); end else if DoTail then Display_time(j=0,Write_time_body, eph_time+JD_MJD,Shift_Time); FirstLast; end; end; if Comet and (TailMod>0) then {producing tail:} begin if (TailWait>0) and (j mod TailWait = 0) then begin debug_ch:=readkey; case UpCase(debug_ch) of 'Q','X',#27,#0: begin gm_close; halt end; 'R': TailWait:=0; { 'S': begin StopTail:=true; exit end; } end end else if keypressed then begin debug_ch:=readkey; case UpCase(debug_ch) of 'Q','X',#27,#0: begin gm_close; halt end; { 'S': begin StopTail:=true; exit end; } end; end; xo[j_ori][1]:=xe; xo[j_ori][2]:=ye; xo[j_ori][3]:=ze; ro[j_ori]:=r; new(p_posi); with p_posi^ do begin next_po:=p_start_po; eph_time_po:=eph_time; ST_po:=ST; end; p_start_po:=p_posi; if j>1 then {particles into tail} begin jori_m1:=succ(succ(j_ori)) mod 3 ; r3:=sqr(ro[jori_m1])*ro[jori_m1]; jori_m2:=succ(j_ori) mod 3; v_ej:=vDust[1]/ro[jori_m1]; {kin. energy prop. to illumination?} fa_ini:=Faint_abs_cmag+Faint_add+ RoundUp(250*(Comet_exponent-2)*ln(ro[jori_m1])/l10); for jcoord:=1 to 3 do begin VVel[jcoord]:=SigStep*(xo[j_ori][jcoord]-xo[jori_m2][jcoord])/2; VVelEj[jcoord]:=-xo[jori_m1][jcoord]*v_ej; VAcc[jcoord]:=-xo[jori_m1][jcoord]*GSol/r3 end; if Rot_angle<>0 then begin MatV(MRot,TailPole,VVelEj); MatVinv(MVel,TailPole,VVelEj); Mat_vect(MRot,VVelEj); Mat(Mrot,3,-Rot_angle); Mat_vect(MRot,VVelEj); Mat_vect(MVel,VVelEj); end; if Tilts then if PoleGiven then MatVinv(MVel,VVelEj,TailPole) else MatVinv(MVel,VVelEj,VTrajPole); for j_tilt:=1 to n_Tilts do for j_dir:=0 to PartTilt[j_tilt].n_dir-1 do begin with PartTilt[j_tilt] do if tiltr>0 then begin P3s_vector(pi/2+2*pi*j_dir/n_dir,pi/2-tiltr,VVelEj); Mat_vect(MVel,VVelEj); for jcoord:= 1 to 3 do VVelEj[jcoord]:=VVelEj[jcoord]*vDust[1]; end; for jjj:=1 to n_Acc do with p_part^ do begin for jcoord:= 1 to 3 do begin xp[jcoord]:=xo[jori_m1][jcoord]+VVel[jcoord] +VVelEj[jcoord]*vDust[jjj]/vDust[1] +VAcc[jcoord]/2; vp[jcoord]:=VVel[jcoord] +VVelEj[jcoord]*vDust[jjj]/vDust[1] +VAcc[jcoord]; end; old:=1; kin:=jjj; central:=(j_tilt=1) and (j_dir=0); fa:=fa_ini; if jjj=n_Acc then dec(fa,100); p_posi:=p_start_po; p_posi:=p_posi^.next_po; p_posi:=p_posi^.next_po; ST:=p_posi^.St_po; eph_time_comp:=p_posi^.Eph_time_po; if (j-2) mod TailMod = 0 then ShowTailParticle; if j > 2 then for j_old:=3 to j+Next do begin p_posi:=p_posi^.next_po; ST:=p_posi^.St_po; eph_time_comp:=p_posi^.Eph_time_po; ShiftOldParticle; end; p_part:=p_start; if p_part^.nextP<> nil then {disposing splitted particles} begin p_part_ins:=p_part^.nextP; p_part^.nextP:=nil; repeat p_part_new:=p_part_ins; p_part_ins:=p_part_new^.nextP; dispose(p_part_new); until p_part_ins=nil; end; end; {of for jjj:=1 to n_Acc do with p_part^ do} end; {of for j_tilt:=1 to n_Tilts do for j_dir:=0 to PartTilt[j_tilt].n_dir-1 do} end; {of particles into tail} j_ori:=succ(j_ori) mod 3; inc(NoPlotted); progress_report; end; {of producing tail} NewTime; end; if DoTail then repeat p_part_new:=p_part; p_part:=p_part^.nextP; dispose(p_part_new); until p_part=nil; if BeforeComputing and (p_start_po<>nil) then begin p_posi:=p_start_po; repeat p_posi_new:=p_posi; p_posi:=p_posi^.next_po; dispose(p_posi_new); until p_posi=nil; end; end; {of NextBefore} begin {of the do_body} if NewMJD then Get_JD else MJD:=JD-JD_MJD; MJD_start:=MJD; if Morning_or_evening<>' ' then begin step:=trunc(step); if step=0 then step:=1; end; if Perih_comp_length<>-999 then perih_from_asc_node:=perih_comp_length-asc_node; if asc_node=-999 then help_p('Longitude of ascending node (an) is missing.'); if perih_from_asc_node=-999 then help_p('Elongation of perihelion from ascending node (pn) is missing.'); if inclin=-999 then help_p('Inclination (i) is missing.'); if semiaxis_large=-999 then if perihel_dist=-999 then help_p('Both large semiaxis (a) and perihelion distance (q) are missing.') else semiaxis_large:=0; if perihel_dist=-999 then perihel_dist:=0; if perih_pass=-999 then if anomaly_date>-998 then if mean_anomaly>-998 then begin if (daily_motion>-998) and (daily_motion<>0) then perih_pass:= anomaly_date-(mean_anomaly/daily_motion) else if semiaxis_large>0 then perih_pass:= anomaly_date- (mean_anomaly*sqrt(semiaxis_large*sqr(semiaxis_large))*365.255/360) else perih_pass:= anomaly_date- (mean_anomaly*sqrt(perihel_dist/(1-excentricity) *sqr(perihel_dist/(1-excentricity)))*365.255/360) end else help_p( 'Time of perihelion passage (pM) (or mean anomaly and its date) is missing.'); if (TailBefore>0) and (TailMod=0) then TailMod:=TailBefore; DoTail:= comet and ((Next>1) or (Before>1) or (TailBefore>1)) and (TailMod>0); { StopTail:=false;} disp_times:= (round(10*step) mod 10 <> 0) or (Morning_or_evening<>' ') or (step<1); eph_time:=MJD; gm_comm('Position(s) of a body according to file '+elm_file_name); if not DoTail then begin N_other:=Next; if N_other>0 then NextBefore; end else begin if Next>0 then begin inc(Before,Next); inc(TailBefore,Next); eph_time:=eph_time+step*Next; Next:=0; if Morning_or_evening<>' ' then begin MoE:=Morning_or_evening; if MoE='B' then MoE:='M'; if not At_height(MoE,-Sun_depth,eph_time) then halt; JD:=eph_time+JD_MJD; Name:=DoW(JD,t_UT,_English); DD:=ND; Get_JD; if surface =cylinder then ST:=SM; end; end; gm_comm('and its tail:'); end; BeforeComputing:=true; N_other:=Before; if TailBefore>N_other then N_other:=TailBefore; if DoTail and Monoscopic then begin moveto(5,report_h); heading:='days/all_before, free mem.: '; outtext(heading); xnum:=getx; NoS:=N_other; NoPlotted:=0; progress_report; end; NextBefore; if DoTail and Monoscopic then begin setcolor(0); moveto(5,report_h); outtext(heading+nst); nst:=''; setcolor(wrtcolour); end; Symb:='+'; NewMJD:=false; SM:=OSM; if (surface=cylinder ) and (Morning_or_evening<>' ') then ST:=SM; JD:=OJD; NoPlotted:=0; end; {of the do_body} function New_naked_Par(j:byte):boolean; begin Par:=ParamStr(j); if Par[1] in ['-','/'] then Par:=copy(Par,2,length(Par)-1); if Par[1] in rem_char then remark:=true; New_naked_Par:= not remark; end; Procedure elm_file(elm_format:char); begin elm_file_Name:=copy(par,3,length(par)-2); Par:=elm_file_name; if elm_file_Name<>old_elm_file_Name then begin inc(Elm_file_No); if Body_No1 then ElmCatNameDefault:=ItemStrD(':',2,Par); if pos('mpc',ElmCatNameDefault)>0 then ElmCatNameKind:=2; if not get_orbital_elements( ElmCatNameKind, ElmCatNameDefault, ItemStrD(':',1,Par)) then help_p(ItemStrD(':',1,Par)+' not found'); Faint_abs_cmag:=roundup(Faint_abs*100); end; end; Found_Body:=true; end; old_elm_file_Name:=elm_file_Name; end; end; Procedure Find_next_elm_file; var par_j,i_pl,jj:byte; begin Elm_file_No:=0; Found_Body:=false; remark:=false; par_j:=1; while (par_j<= paramcount) and (not Found_body) do begin if New_naked_Par(par_j) and (Par[1] = '@') and (length(Par)>1) then begin par_file_name:=copy(par,2,length(par)-1); if not rt_reset_suf(Loc_Obj,par_file,par_file_name,'cnf') then help_p('@ file >>'+par_file_name+'<< not found'); while (not eof(par_file)) and (not Found_body) do begin readln(par_file,par_line); remark_line_file:=false; i_pl:=itemcount(par_line); jj:=1; while (jj<= i_pl) and (not Found_body) and (not remark_line_file) do begin Par:=itemstr(jj,par_line); if Par[1] in rem_char then remark_line_file:=true; if not remark_line_file then if (Par[1]='$') and (length(Par)>2) and (Par[2] in ['!','@']) then elm_file(Par[2]); inc(jj); end; end; close(par_file); end; inc(Par_j); end; if not Found_Body then begin remark:=false; for par_j:= 1 to paramcount do begin if New_naked_Par(par_j) and (Par[1]='$') and (length(Par)>2) and (Par[2] in ['!','@']) then elm_file(Par[2]) end; end; end; Procedure Par_before_follow_elm_file; {in all @files and then in the command-line} var par_j,i_pl,jj:byte; { ContinuePar:boolean;} old_elm_file_Name_loc:Pathstr; const ContinuePar:boolean=true; begin Elm_file_No:=0; old_elm_file_Name_loc:=''; remark:=false; for par_j:= 1 to paramcount do begin if New_naked_Par(par_j) and (Par[1] = '@') and (length(Par)>1) then begin par_file_name:=copy(par,2,length(par)-1); if not rt_reset_suf(Loc_Obj,par_file,par_file_name,'cnf') then help_p('@ file >>'+par_file_name+'<< not found'); ContinuePar:=true; while (not eof(par_file)) and ContinuePar do begin readln(par_file,par_line); remark_line_file:=false; i_pl:=itemcount(par_line); for jj:= 1 to i_pl do begin Par:=itemstr(jj,par_line); if Par[1] in rem_char then remark_line_file:=true; if not remark_line_file then if (Par[1]='$') then if (Par[2] in ['@','!']) then begin elm_file_Name:=copy(par,3,length(par)-2); if elm_file_Name<>old_elm_file_Name_loc then begin inc(Elm_file_No); if Elm_file_No>Body_no then ContinuePar:=false; old_elm_file_Name_loc:=elm_file_Name end end else if ContinuePar then begin Par:=copy(Par,2,length(Par)-1); one_p_param(false); end; end; end; close(par_file); end; end; remark:=false; for par_j:= 1 to paramcount do begin if New_naked_par(par_j) and (Par[1]='$') and (length(Par)>1) then case Par[2] of '@','!': begin inc(Elm_file_No); if Elm_file_No>Body_no then ContinuePar:=false; end; else if ContinuePar then begin Par:=copy(Par,2,length(Par)-1); one_p_param(false); end; end; end; end; begin {of the posi_elm} if paramcount=0 then exit; for i:= 1 to paramcount do if New_naked_par(i) and (Par[1]='$') and (length(Par)>1) then case Par[2] of '?': help_p(''); 'H','h': help_p_file; end; gm_Star_full_mode:=OrPut; Init_next_bef; repeat Find_next_elm_file; {in all @files and then in the command-line, if found, then Read_elm_file } if not Found_Body then exit; Par_before_follow_elm_file; {in all @files and then in the command-line} Do_body; until false; end; {of posi_elm} procedure plan_moon(Jd_p,sm_p:double;n_plp:byte;Moon_only:boolean); var S_L,ra_sol_d,de_sol_d,an,A_noprec,h_noprec:real; j:word; sol_aradp, moon_aradp,Saxe_B:word; DescIsSet:boolean; xmoon,ymoon:longint; namePl:string; emptyPl:boolean; const MSS:word=100; {minimum radius of Sun and Moon in pixels*16} moon_rad:double=0.27248; {in Earth radii} procedure ChangeStereoColorRight; begin if StereoColor then begin if Gray16_stars then begin StereoRight:=true; setcolor(wrtcolour shl TwoOrFour) end else SetColor(right_eye); if GM_mode=3 then gm_verbatim(' pcr '); end; end; procedure ChangeStereoColorLeft; begin if StereoColor then begin if Gray16_stars then begin StereoRight:=false; setcolor(wrtcolour) end else SetColor(lightred); if GM_mode=3 then gm_verbatim(' pcl '); end; end; procedure Par_gm_circle_empty(xc,yc,radius:longint;Lab_st:string); begin if stereoscopic then begin xc:=(xc div 16)*16 + 11; {for getting smooth one-pixel screen parallax step, but coarser stellar positions in PostScript file as well} x_left:=xc+Parallax; if StereoMirror then x_left:= StereoHalfX - x_left; gm_circle_empty(x_left ,yc,radius); if lab_st<>'' then gm_outxy(x_left,yc,lab_st,'c',round_rot(A)); ChangeStereoColorRight; x_left:=xc+stereo_shift-Parallax; gm_circle_empty(x_left ,yc,radius); if lab_st<>'' then gm_outxy(x_left,yc,lab_st,'c',round_rot(A)); ChangeStereoColorLeft; end else begin gm_circle_empty(xc,yc,radius); if lab_st<>'' then gm_outxy(xc,yc,lab_st,'c',round_rot(A)); end; end; procedure Par_gm_circle_full(xc,yc,radius:longint); begin if stereoscopic then begin xc:=(xc div 16)*16 + 11; {for getting smooth one-pixel screen parallax step, but coarser stellar positions in PostScript file as well} x_left:=xc+Parallax; if StereoMirror then x_left:= StereoHalfX - x_left; gm_circle_full(x_left ,yc,radius); ChangeStereoColorRight; gm_circle_full(xc+stereo_shift-Parallax ,yc,radius); ChangeStereoColorLeft; end else gm_circle_full(xc,yc,radius); end; procedure Par_gm_ellipse(xc,yc:longint; anrot,anbeg,anend:double; short,long:word); begin if stereoscopic then begin xc:=(xc div 16)*16 + 11; {for getting smooth one-pixel screen parallax step, but coarser stellar positions in PostScript file as well} x_left:=xc+Parallax; if StereoMirror then x_left:= StereoHalfX - x_left; gm_ellipse(x_left ,yc,anrot,anbeg,anend,short,long); ChangeStereoColorRight; gm_ellipse(xc+stereo_shift-Parallax ,yc,anrot,anbeg,anend,short,long); ChangeStereoColorLeft end else gm_ellipse(xc,yc,anrot,anbeg,anend,short,long); end; procedure Par_gm_line(xl,yl,xr,yr:longint); begin if stereoscopic then begin xl:=(xl div 16)*16 + 11; {for getting smooth one-pixel screen parallax step, but coarser stellar positions in PostScript file as well} x_left:=xl+Parallax; if StereoMirror then x_left:= StereoHalfX - x_left; gm_line(x_left ,yl,xr+Parallax+8,yr); ChangeStereoColorRight; gm_line(xl+stereo_shift-Parallax,yl,xr+stereo_shift-Parallax,yr); ChangeStereoColorLeft; end else gm_line(xl,yl,xr,yr); end; begin DescIsSet:=false; Planet_hel(3,jd_p+tt_ut,RA,Decl,r,name); S_L:=RA+180; Ecl_to_Equ(S_L,ra_d,de_d); if not PrecFixedStars then Precess(ra_d,de_d); ra_sol_d:=ra_d; de_sol_d:=de_d; if not Moon_Only then begin RADe2xy; if aux_boo then begin gm_comm('circle at the Sun position:'); sol_aradp:=round(p16_p_d*rad_to_deg(sol_rad/r)); if sol_aradp3 then begin Planet_hel(j,jd_p+tt_ut,RA,Decl{l,b in fact},r,name); Hel_to_Geo(jd_p+tt_ut,RA,Decl,r,lg,bg,rg); LB_RD(lg,bg,RA_d,De_d); if not PrecFixedStars then Precess(ra_d,de_d); RADe2xy; if aux_boo then begin V:=round(100*planet_faintness (j,object_phase_angleD(RA,Decl,lg,bg),r,rg)); if Elevate_by_refraction then inc(V,round(ZenExt*Air_Mass)); if (V <= Faintness_Lim) or (not PlanFullMode) then begin if V2RV0) then if (n_plp mod WtTimeModMoon) = 0 then begin if ShowDistances then name:=';'+SRx(1,r)+'r.E.' else name:=''; if First_plan[10] or not DescIsSet then begin Display_time(First_plan[10],WtTimeModMoon>0,jd_p,moon_aradp); First_plan[10]:=false; end else gm_outxtxy(x,y,'',Day+Time+name,'','l',round_rot(A), moon_aradp+16,0); end; end; FirstLast end; A:=A_noprec; h:=h_noprec; if h>-0.04 then begin r:=r-sin(h); Parallax:=round(stereo_coef_Moon/r); moon_aradp:=round(p16_p_d*rad_to_deg(moon_rad/r)); if moon_aradp180 then an:=an+180; an:=Bas_Int(-an,360); x:=xmoon; y:=ymoon; {for PS, ellipses are commented out and verbatim arcs are used instead, to have a nice light gray result:} if GM_mode=3 then GM_verbatim('%'); if j<181 then begin {growing moon} Par_gm_ellipse(x,y,an,270,90,moon_aradp,moon_aradp); if GM_mode=3 then gm_verbatim('gsave'+cl+ SI(5,x)+SI(5,MaxY-y)+' translate'+SRe(6,1,an)+ ' rotate'+cl+'0 0 '+SI(5,moon_aradp)+' 270 90 arc'+cl) end else {vaning} begin Par_gm_ellipse(x,y,an,90,270,moon_aradp,moon_aradp); if GM_mode=3 then gm_verbatim('gsave'+cl+ SI(5,x)+SI(5,MaxY-y)+' translate'+SRe(6,1,an)+ ' rotate'+cl+'0 0 '+SI(5,moon_aradp)+' 90 270 arc'+cl) end; if GM_mode=3 then GM_verbatim('%'); case j of 0.. 88, {before first quarter} 181..268: begin {full to last} Par_gm_ellipse(x,y,an,270,90,SAxe_b,moon_aradp); if GM_mode=3 then if j<89 then gm_verbatim(SI(5,SAxe_b)+SI(5,moon_aradp)+' div -1 scale'+cl+ '0 0 '+SI(5,moon_aradp)+' 270 90 arc') else gm_verbatim(SI(5,SAxe_b)+SI(5,moon_aradp)+' div 1 scale'+cl+ '0 0 '+SI(5,moon_aradp)+' 270 90 arc') end; 89.. 91, {first quarter} 269..271: begin {last} Par_gm_line(x,y-moon_aradp,x,y+moon_aradp); if GM_mode=3 then if j<92 then gm_verbatim('0 '+SI(5,-moon_aradp*2)+' rlineto'+cl) else gm_verbatim('0 '+SI(5,moon_aradp*2)+' rlineto'+cl) end; 92..180, {first to full} 272..360: begin {vaning} Par_gm_ellipse(x,y,an,90,270,SAxe_b,moon_aradp); if GM_mode=3 then if j<181 then gm_verbatim(SI(5,SAxe_b)+SI(5,moon_aradp)+' div 1 scale'+cl+ '0 0 '+SI(5,moon_aradp)+' 90 270 arc') else gm_verbatim(SI(5,SAxe_b)+SI(5,moon_aradp)+' div -1 scale'+cl+ '0 0 '+SI(5,moon_aradp)+' 90 270 arc') end; end; if GM_mode=3 then gm_verbatim(cl+ ' gsave 0.85 setgray fill grestore'+cl+ ' stroke'+cl+ 'grestore'+cl); if not StereoColor then setcolor(wrtcolour); end; {of Moon over horizon} end; {of the Moon part} end; {of plan_moon} end.