c shift to Linux 25/4/05: c rm'ed ifix line 1344 c f77 -fno-backslash -o prc prc.f -lpgplot -L/usr/X11R6/lib -lX11 c c PRC v1.0 5/96 c v1.0.1 1/6/96 bug in Mach(beta) fixed c v1.0.2 1/4/97 bug in backspace fixed c v1.0.3 1/8/97 lenj function made portable c NB. Reset '/xw' in pgbegin, modify alternate form in mread, as needed. c Under SunOS/Solaris: f77 -o prc -xl prc.f -lpgplot -lX11 c NB. DOS USERS: contact author if you are interested in a DOS version c that uses the true mouse, rather than pgplot's arrow keys. c c Usage: "click" means use the mouse (Suns etc.) or arrow keys (PC) c to move the cursor, and click a mouse button or press any key. c o click on a function button; c o type the required value in each entry box as prompted, finishing c each with the ENTER/RETURN key. c o Note -- each field retains its previous value if any; ENTER/RETURN c to use, or DEL/BACKSPACE or click on box label to clear. c o Note -- as with all dialog box-based utilities, the cursor must c be in the box to enter values. c c NB. A brief description of button functions read left to right, bottom c to top, precedes their subroutines. c real mu,kaxis character*35 cbut(36) character*14 stin1,stin2,stin3,stin4,stout1,stout2,stout3 character*14 labin1,labin2,labin3,labin4,labout1,labout2,labout3 character*26 qis,qask,qid c nb. overide to vpright in pggett/pgputt c nb. force curs location in pggett c nb. sleep off qis='PRC v. 1.0.2' qask='Comments/Questions:' qid='hughes@astro.lsa.umich.edu' c call pgbegin (0,'/xw',1,1) call pgask (.false.) c cbut(01)='\gg(\gb)' cbut(02)='\gb(\gg)' cbut(03)='D(\gg,\gh\do\u)' cbut(04)='D(\gb,\gh\do\u)' cbut(05)='\gb\da\u(\gg,\gh\do\u)' cbut(06)='\gb\da\u(\gb,\gh\do\u)' cbut(07)='\gh\do\u(\gh\di\u,\gg)' cbut(08)='\gh\do\u(\gh\di\u,\gb)' cbut(09)='\gh\di\u(\gh\do\u,\gg)' cbut(10)='\gh\di\u(\gh\do\u,\gb)' cbut(11)='\gh\do\u(D,\gg)' cbut(12)='\gh\do\u(D,\gb)' cbut(13)='\gh\do\u(\gb\da\u,\gg)' cbut(14)='\gh\do\u(\gb\da\u,\gb)' cbut(15)='\gg(D,\gh\do\u)' cbut(16)='\gb(D,\gh\do\u)' cbut(17)='\gg(\gb\da\u,\gh\do\u)' cbut(18)='\gb(\gb\da\u,\gh\do\u)' cbut(19)='\gb,\gg,\gh\do\u(D,\gb\da\u)' cbut(20)='\gb\da\u\umax\d,\gh\do\u,D(\gg)' cbut(21)='\gb\da\u\umax\d,\gh\do\u,D(\gb)' cbut(22)='R\d\gh\u(z,q0)' cbut(23)='R\dl\u(z,q0)' cbut(24)='\gb\da\u(z,q0,H0,\gm)' cbut(25)='\gm(z,q0,H0,\gb\da\u)' cbut(26)='D\dmax\u(\gb)' cbut(27)='D\dmax\u(\gg)' cbut(28)='P(\gd,\gq,k)' cbut(29)='k(\gd,\gq,P)' cbut(30)='\gq(\gd,k,P)' cbut(31)='\gb+\gb' cbut(32)='\gg+\gg' cbut(33)='M(\gg,\gg\ds\u)' cbut(34)='M(\gb,\gb\ds\u)' cbut(35)=' ?' cbut(36)='Quit' c padl=0.05 padr=0.95 padd=0.05 padh=0.6 vpinleft=0.05 vpinhite=0.0615 vpin1=1.0-vpinhite vpin2=vpin1-1.5*vpinhite vpin3=vpin2-1.5*vpinhite vpin4=vpin3-1.5*vpinhite vpoutleft=0.55 vpouthite=0.0615 vpout1=1.0-vpouthite vpout2=vpout1-1.5*vpouthite vpout3=vpout2-1.5*vpouthite vpout4=vpout3-1.5*vpouthite c vpinret1=0.5 vpinret2=0.5 vpinret3=0.5 vpinret4=0.5 vpoutret1=0.9 vpoutret2=0.9 vpoutret3=0.9 vpoutret4=0.9 c index=0 call pgbut (padl,padd,padr,padh,4,9,cbut,index) if (index.lt.0) then write (*,*) 'oops!' stop end if index=1 !arb. value 20 call pgbut (padl,padd,padr,padh,4,9,cbut,index) if (index.lt.0) then write (*,*) 'oops!' stop end if call pgdiax (vpinleft,vpin1,vpinhite,vpinret1) call pgdiax (vpinleft,vpin2,vpinhite,vpinret2) call pgdiax (vpinleft,vpin3,vpinhite,vpinret3) call pgdiax (vpinleft,vpin4,vpinhite,vpinret4) go to (30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180, *190,200,210,220,230,240,250,260,320,270,280,290,300,310,330,340, *350,360,370,380),index c 30 labin1='\gb:' 31 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*31) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(beta).ge.1.0) then stout1='invalid' else x=g(beta) write (stout1,'(f10.5)') x end if labout1='\gg:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 40 labin1='\gg:' 41 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*41) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0) then stout1='invalid' else x=b(gamma) write (stout1,'(f10.5)') x end if labout1='\gb:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 50 labin1='\gg:' labin2='\gh:' 51 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*51) 52 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*52) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=dg(gamma,t) write (stout1,'(f10.5)') x end if labout1='D:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 60 labin1='\gb:' labin2='\gh:' 61 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*61) 62 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*62) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (beta.ge.1.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=db(beta,t) write (stout1,'(f10.5)') x end if labout1='D:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 70 labin1='\gg:' labin2='\gh:' 71 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*71) 72 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*72) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=bappg(gamma,t) write (stout1,'(f10.5)') x end if labout1='\gb\da\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 80 labin1='\gb:' labin2='\gh:' 81 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*81) 82 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*82) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (beta.ge.1.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=bappb(beta,t) write (stout1,'(f10.5)') x end if labout1='\gb\da\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 90 labin1='\gh:' labin2='\gg:' 91 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,ti,*91) 92 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,gamma,*92) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (ti.lt.0.0.or.ti.gt.180.0.or.gamma.lt.1.0) then stout1='invalid' else x=tg(ti,gamma) write (stout1,'(f10.5)') x end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 100 labin1='\gh:' labin2='\gb:' 101 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,ti,*101) 102 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,beta,*102) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (ti.lt.0.0.or.ti.gt.180.0.or.abs(beta).ge.1.0) then stout1='invalid' else x=tb(ti,beta) write (stout1,'(f10.5)') x end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 110 labin1='\gh:' labin2='\gg:' 111 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,ti,*111) 112 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,gamma,*112) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (ti.lt.0.0.or.ti.gt.180.0.or.gamma.lt.1.0) then stout1='invalid' else x=tig(ti,gamma) write (stout1,'(f10.5)') x end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 120 labin1='\gh:' labin2='\gb:' 121 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,ti,*121) 122 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,beta,*122) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (ti.lt.0.0.or.ti.gt.180.0.or.abs(beta).ge.1.0) then stout1='invalid' else x=tib(ti,beta) write (stout1,'(f10.5)') x end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 130 labin1='D:' labin2='\gg:' 131 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,d,*131) 132 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,gamma,*132) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (d.lt.0.0.or.gamma.lt.1.0) then stout1='invalid' else x=tdg(d,gamma) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 140 labin1='D:' labin2='\gb:' 141 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,d,*141) 142 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,beta,*142) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (d.lt.0.0.or.abs(beta).ge.1.0) then stout1='invalid' else x=tdb(d,beta) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 150 labin1='\gb\da\u:' labin2='\gg:' 151 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,bapp,*151) 152 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,gamma,*152) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0) then stout1='invalid' t1=-0.1e10 else call tbappg (bapp,gamma,t1,t2) if (t1.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') t1 write (stout2,'(f10.5)') t2 end if end if labout1='\gh\d1\u:' labout2='\gh\d2\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (t1.gt.-0.1e1) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) go to 20 160 labin1='\gb\da\u:' labin2='\gb:' 161 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,bapp,*161) 162 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,beta,*162) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(beta).ge.1.0) then stout1='invalid' t1=-0.1e10 else call tbappb (bapp,beta,t1,t2) if (t1.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') t1 write (stout2,'(f10.5)') t2 end if end if labout1='\gh\d1\u:' labout2='\gh\d2\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (t1.gt.-0.1e1) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) go to 20 170 labin1='D:' labin2='\gh:' 171 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,d,*171) 172 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*172) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (d.lt.0.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' g1=-0.1e10 else call gdt (d,t,g1,g2) if (g1.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') g1 write (stout2,'(f10.5)') g2 end if end if labout1='\gg\d1\u:' labout2='\gg\d2\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (g1.gt.-0.1e1) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) go to 20 180 labin1='D:' labin2='\gh:' 181 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,d,*181) 182 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*182) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (d.lt.0.0.or.t.lt.0.0.or.t.gt.180.0) then stout1='invalid' b1=-0.1e10 else call bdt (d,t,b1,b2) if (b1.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') b1 write (stout2,'(f10.5)') b2 end if end if labout1='\gb\d1\u:' labout2='\gb\d2\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (b1.gt.-0.1e1) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) go to 20 190 labin1='\gb\da\u:' labin2='\gh:' 191 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,bapp,*191) 192 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*192) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=gbappt(bapp,t) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='\gg:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 200 labin1='\gb\da\u:' labin2='\gh:' 201 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,bapp,*201) 202 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*202) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (t.lt.0.0.or.t.gt.180.0) then stout1='invalid' else x=bbappt(bapp,t) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='\gb:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 210 labin1='D:' labin2='\gb\da\u:' 211 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,d,*211) 212 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,bapp,*212) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (t.lt.0.0.or.t.gt.180.0) then stout1='invalid' beta=-0.1e10 else call dandb (d,bapp,beta,gamma,t) if (beta.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') beta write (stout2,'(f10.5)') gamma write (stout3,'(f10.5)') t end if end if labout1='\gb:' labout2='\gg:' labout3='\gh:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (beta.gt.-0.1e1) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) if (beta.gt.-0.1e1) *vpoutret3=pgputt (vpoutleft,vpout3,vpouthite,labout3,stout3) go to 20 220 labin1='\gg:' 221 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*221) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0) then stout1='invalid' else call bmaxg (gamma,bapmx,t,d) write (stout1,'(f10.5)') bapmx write (stout2,'(f10.5)') t write (stout3,'(f10.5)') d end if labout1='\gb\da\u\u*\d:' labout2='\gh:' labout3='D:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (gamma.ge.1.0) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) if (gamma.ge.1.0) *vpoutret3=pgputt (vpoutleft,vpout3,vpouthite,labout3,stout3) go to 20 230 labin1='\gb:' 231 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*231) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(beta).ge.1.0) then stout1='invalid' else call bmaxb (beta,bapmx,t,d) write (stout1,'(f10.5)') bapmx write (stout2,'(f10.5)') t write (stout3,'(f10.5)') d end if labout1='\gb\da\u\u*\d:' labout2='\gh:' labout3='D:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) if (abs(beta).lt.1.0) *vpoutret2=pgputt (vpoutleft,vpout2,vpouthite,labout2,stout2) if (abs(beta).lt.1.0) *vpoutret3=pgputt (vpoutleft,vpout3,vpouthite,labout3,stout3) go to 20 240 labin1='z:' labin2='q0:' 241 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,z,*241) 242 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,q0,*242) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (z.lt.0.0.or.q0.lt.0.0) then stout1='invalid' else x=ra(z,q0) write (stout1,'(f10.5)') x end if labout1='R\d\gh\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 250 labin1='z:' labin2='q0:' 251 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,z,*251) 252 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,q0,*252) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (z.lt.0.0.or.q0.lt.0.0) then stout1='invalid' else x=rl(z,q0) write (stout1,'(f10.5)') x end if labout1='R\dl\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 260 labin1='z:' labin2='q0:' labin3='H0:' labin4='\gm:' 261 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,z,*261) 262 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,q0,*262) 263 vpinret3=pggett (vpinleft,vpin3,vpinhite,labin3,stin3) call mread (stin3,h0,*263) 264 vpinret4=pggett (vpinleft,vpin4,vpinhite,labin4,stin4) call mread (stin4,mu,*264) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (z.lt.0.0.or.q0.lt.0.0.or.h0.lt.0.0) then stout1='invalid' else x=bappmu(z,q0,h0,mu) write (stout1,'(f10.5)') x end if labout1='\gb\da\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 270 labin1='\gb:' 271 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*271) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(beta).ge.1.0) then stout1='invalid' else x=dmaxb(beta) write (stout1,'(f10.5)') x end if labout1='D\u*\d:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 280 labin1='\gg:' 281 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*281) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0) then stout1='invalid' else x=dmaxg(gamma) write (stout1,'(f10.5)') x end if labout1='D\u*\d:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 290 labin1='\gd:' labin2='\gq:' labin3='k:' 291 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,delta,*291) 292 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*292) 293 vpinret3=pggett (vpinleft,vpin3,vpinhite,labin3,stin3) call mread (stin3,kaxis,*293) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (delta.le.-1.0.or.t.lt.0.0.or.t.gt.90.0.or.kaxis.lt.0.0.or. *kaxis.gt.1.0) then stout1='invalid' else x=pol(delta,t,kaxis) write (stout1,'(f10.5)') x end if labout1='P:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 300 labin1='\gd:' labin2='\gq:' labin3='P:' 301 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,delta,*301) 302 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,t,*302) 303 vpinret3=pggett (vpinleft,vpin3,vpinhite,labin3,stin3) call mread (stin3,xpol,*303) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (delta.le.-1.0.or.t.lt.0.0.or.t.gt.90.0.or.xpol.lt.0.0.or.xpol *.gt.1.0) then stout1='invalid' else x=axisk(delta,t,xpol) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='k:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 310 labin1='\gd:' labin2='k:' labin3='P:' 311 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,delta,*311) 312 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,kaxis,*312) 313 vpinret3=pggett (vpinleft,vpin3,vpinhite,labin3,stin3) call mread (stin3,xpol,*313) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (delta.le.-1.0.or.xpol.lt.0.0.or.xpol.gt.1.0.or.kaxis.lt.0.0. *or.kaxis.ge.1.0) then stout1='invalid' else x=tpol(delta,kaxis,xpol) if (x.lt.-0.1e5) then stout1='no solution' else write (stout1,'(f10.5)') x end if end if labout1='\gq:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 320 labin1='z:' labin2='q0:' labin3='H0:' labin4='\gb\da\u:' 321 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,z,*321) 322 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,q0,*322) 323 vpinret3=pggett (vpinleft,vpin3,vpinhite,labin3,stin3) call mread (stin3,h0,*323) 324 vpinret4=pggett (vpinleft,vpin4,vpinhite,labin4,stin4) call mread (stin4,bapp,*324) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (z.lt.0.0.or.q0.lt.0.0.or.h0.lt.0.0) then stout1='invalid' else x=umbapp(z,q0,h0,bapp) write (stout1,'(f10.5)') x end if labout1='\gm:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 330 labin1='\gb\d1\u:' labin2='\gb\d2\u:' 331 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,b1,*331) 332 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,b2,*332) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(b1).ge.1.0.or.abs(b2).ge.1.0) then stout1='invalid' else x=bbb(b1,b2) write (stout1,'(f10.5)') x end if labout1='\gb\d3\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 340 labin1='\gg\d1\u:' labin2='\gg\d2\u:' 341 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,g1,*341) 342 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,g2,*342) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (abs(g1).lt.1.0.or.abs(g2).lt.1.0) then stout1='invalid' else x=ggg(g1,g2) write (stout1,'(f10.5)') x end if labout1='\gg\d3\u:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 350 labin1='\gg:' labin2='\gg\ds\u:' 351 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,gamma,*351) 352 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,gammas,*352) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (gamma.lt.1.0.or.gammas.lt.1.0) then stout1='invalid' else x=rmachg(gamma,gammas) write (stout1,'(f10.5)') x end if labout1='M:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 360 labin1='\gb:' labin2='\gb\ds\u:' 361 vpinret1=pggett (vpinleft,vpin1,vpinhite,labin1,stin1) call mread (stin1,beta,*361) 362 vpinret2=pggett (vpinleft,vpin2,vpinhite,labin2,stin2) call mread (stin2,betas,*362) call pgdiax (vpoutleft,vpout1,vpouthite,vpoutret1) call pgdiax (vpoutleft,vpout2,vpouthite,vpoutret2) call pgdiax (vpoutleft,vpout3,vpouthite,vpoutret3) call pgdiax (vpoutleft,vpout4,vpouthite,vpoutret4) if (beta.ge.1.0.or.betas.ge.1.0) then stout1='invalid' else x=rmachb(beta,betas) write (stout1,'(f10.5)') x end if labout1='M:' vpoutret1=pgputt (vpoutleft,vpout1,vpouthite,labout1,stout1) go to 20 370 vpinret1=pgputt (vpinleft,vpin1,vpinhite,' ',qis) vpinret2=pgputt (vpinleft,vpin2,vpinhite,' ',qask) vpinret3=pgputt (vpinleft,vpin3,vpinhite,' ',qid) go to 20 c 380 call pgend () stop end c% c THE FOLLOWING FORTRAN CALLABLE FUNCTIONS PROVIDE QUICK ACCESS c TO THE BASIC COMPUTATIONS OF SPECIAL RELATIVITY USED IN THE c STUDY OF RELATIVISTIC JETS. c c V 1.0 P. A. HUGHES c V 1.1 c TESTS INSERTED INTO TDG, TDB, GBAPPT AND BBAPPT. c V 2.0 c UMBAPP, DMAX, POL AND TRANSFORMATION ROUTINES ADDED. c c NB. APPEND THE ENTIRE PACKET OF FUNCTIONS: SOME ARE c INTERDEPENDENT. c c c beta beta c <-- < -------------:------------- > --> c 0 \ . 0 c \ . c \ . 0 <= theta <= 180 c \< c \ c TO OBSERVER c c NB. INPUT/OUTPUT IS IN DEGREES. c c CONTENTS c ======== c c 1 G(B).........GAMMA AS A FUNCTION OF BETA. c 2 B(G).........BETA AS A FUNCTION OF GAMMA. c 3 DG(G,T)......DOPPLER FACTOR AS A FUNCTION OF GAMMA AND VIEWING ANGLE. c 4 DB(B,T)......DOPPLER FACTOR AS A FUNCTION OF BETA AND VIEWING ANGLE. c 5 BAPPG(G,T)...APPARENT BETA AS A FUNCTION OF GAMMA AND VIEWING ANGLE. c 6 BAPPB(B,T)...APPARENT BETA AS A FUNCTION OF BETA AND VIEWING ANGLE. c 7 TG(TI,G).....VIEWING ANGLE AS A FUNCTION OF INTRINSIC EMISSION ANGLE c AND GAMMA. c 8 TB(TI,B).....VIEWING ANGLE AS A FUNCTION OF INTRINSIC EMISSION ANGLE c AND BETA. c 9 TIG(TI,G)....INTRINSIC EMISSION ANGLE AS A FUNCTION OF VIEWING ANGLE c AND GAMMA. c 10 TIB(TI,B)....INTRINSIC EMISSION ANGLE AS A FUNCTION OF VIEWING ANGLE c AND BETA. c 11 TDG(D,G).....VIEWING ANGLE AS A FUNCTION OF DOPPLER FACTOR AND GAMMA; c -INFINITY IF NO SOLUTION. c 12 TDB(D,B).....VIEWING ANGLE AS A FUNCTION OF DOPPLER FACTOR AND BETA; c -INFINITY IF NO SOLUTION. c 13 TBAPPG(BAPP,G,T1,T2)...SUBROUTINE THAT RETURNS IN T1, T2 THE TWO ANGLES c AT WHICH BAPP IS GIVEN BY GAMMA; -INFINITY IF NO SOLUTION. c 14 TBAPPB(BAPP,B,T1,T2)...SUBROUTINE THAT RETURNS IN T1, T2 THE TWO ANGLES c AT WHICH BAPP IS GIVEN BY BETA; -INFINITY IF NO SOLUTION. c 15 GDT(D,T,G1,G2).........SUBROUTINE THAT RETURNS IN G1, G2 THE TWO GAMMAS c AT WHICH DOPPLER FACTOR IS GIVEN BY ANGLE; -INFINITY IF NO c SOLUTION. c 16 BDT(D,T,B1,B2).........SUBROUTINE THAT RETURNS IN B1, B2 THE TWO BETAS c AT WHICH DOPPLER FACTOR IS GIVEN BY ANGLE; -INFINITY IF NO c SOLUTION. c 17 GBAPPT(BAPP,T)...GAMMA AS A FUNCTION OF BAPP AND VIEWING ANGLE; c -INFINITY IF NO SOLUTION. c 18 BBAPPT(BAPP,T)...BETA AS A FUNCTION OF BAPP AND VIEWING ANGLE; c -INFINITY IF NO SOLUTION. c 19 DANDB(D,BAPP,BETA,G,T)..SUBROUTINE THAT TAKES DOPPLER FACTOR, BAPP AND c RETURNS BETA, GAMMA AND VIEWING ANGLE THAT YIELD THESE; c -INFINITY IF NO SOLUTION. c 20 BMAXG(GAMMA,BAPMX,T,D)...SUBROUTINE THAT FOR GIVEN GAMMA RETURNS MAXIMUM c BAPP, ANGLE AT WHICH THIS OCCURS AND CORRESPONDING DOPPLER. c 21 BMAXB(BETA,BAPMX,T,D)....SUBROUTINE THAT FOR GIVEN BETA RETURNS MAXIMUM c BAPP, ANGLE AT WHICH THIS OCCURS AND CORRESPONDING DOPPLER. c 22 RA(Z,Q0).....COSMOLOGICAL (ANGULAR DIAMETER) DISTANCE AS A FUNCTION OF c REDSHIFT AND DECELERATION PARAMETER, SCALED WITH C/H0 c (I.E. DISTANCE = (C/H0) * RA(Z,Q0)). c 23 RL(Z,Q0).....COSMOLOGICAL (LUMINOSITY) DISTANCE AS A FUNCTION OF c REDSHIFT AND DECELERATION PARAMETER, SCALED WITH C/H0 c (I.E. DISTANCE = (C/H0) * RL(Z,Q0)). c 24 BAPPMU(Z,Q0,H0,MU)...APPARENT BETA COMPUTED FROM A MEASURED PROPER MOTION c OF MU M.A.S./YEAR, WITH H0 IN KM/SEC/MPC. c 25 UMBAPP(Z,Q0,H0,BAPP)...INVERSE OF ABOVE: PROPER MOTION FROM APPARENT BETA. c 26 DMAXB(BETA)...MAXIMUM DOPPLER FACTOR FOR A GIVEN BETA. c 27 DMAXG(GAMMA)...MAXIMUM DOPPLER FACTOR FOR A GIVEN GAMMA. c 28 POL(DELTA,T,KAXIS)...FRACTIONAL POLARIZATION AS A FUNCTION OF DELTA, c ANGLE FROM PLANE OF COMPRESSION AND COMP. FACTOR. c 29 AXISK(DELTA,T,POL)...COMPRESSION THAT GIVES POL FOR GIVEN DELTA AND ANGLE; c -INFINITY IF NO SOLUTION. c 30 TPOL(DELTA,KAXIS,POL)...ANGLE THAT GIVES POL FOR GIVEN DELTA AND c COMPRESSION; -INFINITY IF NO SOLUTION. c 31 BBB(B1,B2)...LORENTZ TRANSFORMATION: B1 (+) B2 = B3 c 32 GGG(G1,G2)...LORENTZ TRANSFORMATION: G1 (+) G2 = G3; A SIGN IS c ASSOCIATED WITH INPUT AND OUTPUT G. c 33 RMACHG(G,GS)..RELATIVISTIC MACH NUMBER AS A FUNCTION OF GAMMA-FLOW AND c GAMMA-SOUND. c 34 RMACHB(B,BS)..RELATIVISTIC MACH NUMBER AS A FUNCTION OF BETA-FLOW AND c BETA-SOUND. c function g (b) g=1.0/sqrt(1.0-b**2) return end function b (g) if (g.le.1.0) go to 10 b=sqrt(1.0-1.0/g**2) go to 20 10 b=0.0 20 return end function dg (g,t) pi=3.14159265 rt=pi*t/180.0 beta=b(g) dg=1.0/(g*(1.0-beta*cos(rt))) return end function db (b,t) pi=3.14159265 rt=pi*t/180.0 gamma=g(b) db=1.0/(gamma*(1.0-b*cos(rt))) return end function bappg (g,t) pi=3.14159265 rt=pi*t/180.0 beta=b(g) bappg=beta*sin(rt)/(1.0-beta*cos(rt)) return end function bappb (b,t) pi=3.14159265 rt=pi*t/180.0 bappb=b*sin(rt)/(1.0-b*cos(rt)) return end function tg (ti,g) pi=3.14159265 rti=pi*ti/180.0 beta=b(g) rtgc=(beta+cos(rti))/(1.0+beta*cos(rti)) rtg=acos(rtgc) tg=180.0*rtg/pi return end function tb (ti,b) pi=3.14159265 rti=pi*ti/180.0 rtbc=(b+cos(rti))/(1.0+b*cos(rti)) rtb=acos(rtbc) tb=180.0*rtb/pi return end function tig (t,g) pi=3.14159265 rt=pi*t/180.0 beta=b(g) rtigc=(cos(rt)-beta)/(1.0-beta*cos(rt)) rtig=acos(rtigc) tig=180.0*rtig/pi return end function tib (t,b) pi=3.14159265 rt=pi*t/180.0 rtibc=(cos(rt)-b)/(1.0-b*cos(rt)) rtib=acos(rtibc) tib=180.0*rtib/pi return end function tdg (d,g) pi=3.14159265 beta=b(g) ctdg=(1.0-1.0/(d*g))/beta if (abs(ctdg).gt.1.0) go to 10 tdgr=acos(ctdg) tdg=180.0*tdgr/pi go to 20 10 tdg=-.1e10 20 return end function tdb (d,b) pi=3.14159265 gamma=g(b) ctdb=(1.0-1.0/(d*gamma))/b if (abs(ctdb).gt.1.0) go to 10 tdbr=acos(ctdb) tdb=180.0*tdbr/pi go to 20 10 tdb=-.1e10 20 return end subroutine tbappg (bapp,g,t1,t2) pi=3.14159265 beta=b(g) test=beta**2*(1.0+bapp**2)-bapp**2 if (test.lt.0.0) go to 10 ct1=(bapp**2+sqrt(test))/(beta*(1.0+bapp**2)) ct2=(bapp**2-sqrt(test))/(beta*(1.0+bapp**2)) t1r=acos(ct1) t2r=acos(ct2) t1=180.0*t1r/pi t2=180.0*t2r/pi go to 20 10 t1=-.1e10 t2=t1 20 return end subroutine tbappb (bapp,b,t1,t2) pi=3.14159265 beta=b test=beta**2*(1.0+bapp**2)-bapp**2 if (test.lt.0.0) go to 10 ct1=(bapp**2+sqrt(test))/(beta*(1.0+bapp**2)) ct2=(bapp**2-sqrt(test))/(beta*(1.0+bapp**2)) t1r=acos(ct1) t2r=acos(ct2) t1=180.0*t1r/pi t2=180.0*t2r/pi go to 20 10 t1=-.1e10 t2=t1 20 return end subroutine gdt (d,t,g1,g2) pi=3.14159265 tr=t*pi/180.0 test=1.0-d**2*(sin(tr))**2 if (test.lt.0.0) go to 10 b1=(d**2*cos(tr)+sqrt(test))/(1.0+(d*cos(tr))**2) b2=(d**2*cos(tr)-sqrt(test))/(1.0+(d*cos(tr))**2) g1=g(b1) g2=g(b2) go to 20 10 g1=-.1e10 g2=g1 20 return end subroutine bdt (d,t,b1,b2) pi=3.14159265 tr=t*pi/180.0 test=1.0-d**2*(sin(tr))**2 if (test.lt.0.0) go to 10 b1=(d**2*cos(tr)+sqrt(test))/(1.0+(d*cos(tr))**2) b2=(d**2*cos(tr)-sqrt(test))/(1.0+(d*cos(tr))**2) go to 20 10 b1=-.1e10 b2=b1 20 return end function gbappt (bapp,t) pi=3.14159265 tr=t*pi/180.0 beta=bapp/(sin(tr)+bapp*cos(tr)) if (beta.gt.1.0) go to 10 gbappt=g(beta) go to 20 10 gbappt=-.1e10 20 return end function bbappt (bapp,t) pi=3.14159265 tr=t*pi/180.0 bbappt=bapp/(sin(tr)+bapp*cos(tr)) if (bbappt.gt.1.0) bbappt=-.1e10 return end subroutine dandb (d,bapp,beta,g,t) pi=3.14159265 g=d*(1.0+(1.0+bapp**2)/d**2)/2.0 if (g.lt.1.0) go to 10 beta=b(g) ctr=(1.0-1.0/(d*g))/beta if (abs(ctr).gt.1.0) go to 10 tr=acos(ctr) t=180.0*tr/pi go to 20 10 beta=-.1e10 g=beta t=beta 20 return end subroutine bmaxg (gamma,bapmx,t,d) pi=3.14159265 beta=b(gamma) bapmx=beta*gamma tr=acos(beta) t=180.0*tr/pi d=dg(gamma,t) return end subroutine bmaxb (beta,bapmx,t,d) pi=3.14159265 gamma=g(beta) bapmx=gamma*beta tr=acos(beta) t=180.0*tr/pi d=dg(gamma,t) return end function ra (z,q0) ohm=2.0*q0 fact=sqrt(1.0+ohm*z) ratio=(1.0+fact+z)/(1.0+fact+ohm*z/2.0) ra=z*ratio/(1.0+z)**2 return end function rl (z,q0) dist=ra(z,q0) rl=dist*(1.0+z)**2 return end function bappmu (z,q0,h0,mu) real*4 mu dist=ra(z,q0) h=h0/100.0 bappmu=47.4*dist*mu*(1.0+z)/h return end function umbapp (z,q0,h0,bapp) dist=ra(z,q0) h=h0/100.0 umbapp=bapp*h/(47.4*dist*(1.0+z)) return end function dmaxb (beta) dmaxb=sqrt((1.0+beta)/(1.0-beta)) return end function dmaxg (gamma) beta=b(gamma) dmaxg=sqrt((1.0+beta)/(1.0-beta)) return end function pol (delta,t,kaxis) real*4 kaxis pi=3.14159265 rt=pi*t/180.0 f=(delta+1.0)/(delta+7.0/3.0) temp1=1.0-kaxis*kaxis temp2=cos(rt)*cos(rt) pol=f*temp1*temp2/(2.0-temp1*temp2) return end function axisk (delta,t,pol) pi=3.14159265 rt=pi*t/180.0 f=(delta+1.0)/(delta+7.0/3.0) temp=1.0-2.0*pol/((pol+f)*cos(rt)*cos(rt)) if (temp.lt.0.0) go to 10 axisk=sqrt(temp) go to 20 10 axisk=-.1e10 20 return end function tpol (delta,kaxis,pol) real*4 kaxis pi=3.14159265 f=(delta+1.0)/(delta+7.0/3.0) temp=sqrt(2.0*pol/((pol+f)*(1.0-kaxis*kaxis))) if (temp.gt.1.0) go to 10 tpol=acos(temp) tpol=180.0*tpol/pi go to 20 10 tpol=-.1e10 20 return end function bbb (b1,b2) bbb=(b1+b2)/(1.0+b1*b2) return end function ggg (g1,g2) b1=sign(b(abs(g1)),g1) b2=sign(b(abs(g2)),g2) b3=(b1+b2)/(1.0+b1*b2) ggg=sign(g(b3),b3) return end function rmachg (g,gs) beta=b(g) betas=b(gs) rmachg=g*beta/(gs*betas) return end function rmachb (b,bs) gamma=g(b) gammas=g(bs) rmachb=gamma*b/(gammas*bs) return end c% subroutine mread (string,value,*) character*1 ding character*(*) string c---- list directed internal read is not in general portable read (string,*,err=10) value c---- c open (unit=1,status='scratch') c write (1,*) string c rewind (1) c read (1,*,iostat=ibad,err=20) value c20 close (1) c if (ibad.gt.0) go to 10 c---- return 10 string=' ' ding=char(07) write (*,*) ding return 1 end c% subroutine pgbut (padllx,padlly,padurx,padury,n,m,cs,io) real padllx,padlly,padurx,padury ! BB of button pad integer n,m ! nxm array of buttons character*(*) cs(n*m) ! nxm labels for buttons integer io ! 0 to set, >0 to recall integer pgcurs character*1 ccurs real xbox(4),ybox(4) save xcurs,ycurs c if (n.le.0.or.m.le.0.or.n.gt.99.or.m.gt.99) then io=-1 write (*,*) 'PGBUT: bad button array size' return end if if (padllx.lt.0.0.or.padurx.gt.1.0.or.padlly.lt.0.0.or.padury.gt. $ 1.0) then io=-1 write (*,*) 'PGBUT: bad viewport values' return end if if (io.lt.0) then io=-1 write (*,*) 'PGBUT: bad IO value' return end if c call pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c chargoal=0.55 ! in world coords surround=0.45 ! in world coords xtoffset=0.075 ! in world coords ytoffset=0.3 ! in world coords c call pgsvp (padllx,padurx,padlly,padury) call pgswin (0.0,float(n),0.0,float(m)) call pgscf (2) call pgsls (1) call pgslw (2) call pgsch (1.0) call pgqtxt (1.0,1.0,0.0,0.0,'DUMMY',xbox,ybox) boxbase=amin1(ybox(1),ybox(2),ybox(3),ybox(4)) boxtop =amax1(ybox(1),ybox(2),ybox(3),ybox(4)) charhite=boxtop-boxbase charfact=chargoal/charhite call pgsch (charfact) c if (io.eq.0) then call pgsci (0) call pgsfs (1) call pgrect (0.0,float(n),0.0,float(m)) call pgsci (1) call pgsfs (2) call pgrect (0.0,float(n),0.0,float(m)) do i=1,n do j=1,m call pgrect (float(i)-0.5-surround,float(i)-0.5+surround, $ float(j)-0.5-surround,float(j)-0.5+surround) call pgptxt (float(i-1)+xtoffset,float(j-1)+ytoffset,0.0, $ 0.0,cs((j-1)*n+i)) end do end do else j=((io-1)/n)+1 i=io-n*(j-1) call pgsci (0) call pgsfs (1) call pgrect (float(i)-0.5-surround,float(i)-0.5+surround, $ float(j)-0.5-surround,float(j)-0.5+surround) call pgsci (1) call pgsfs (2) call pgrect (float(i)-0.5-surround,float(i)-0.5+surround, $ float(j)-0.5-surround,float(j)-0.5+surround) call pgptxt (float(i-1)+xtoffset,float(j-1)+ytoffset,0.0,0.0 $ ,cs((j-1)*n+i)) ido=pgcurs (xcurs,ycurs,ccurs) i=ifix(xcurs)+1 j=ifix(ycurs)+1 io=(j-1)*n+i call pgsci (4) call pgsfs (1) call pgrect (float(i)-0.5-surround,float(i)-0.5+surround, $ float(j)-0.5-surround,float(j)-0.5+surround) call pgstbg (4) call pgsci (1) call pgptxt (float(i-1)+xtoffset,float(j-1)+ytoffset,0.0,0.0 $ ,cs((j-1)*n+i)) c call sleep (1) end if c call pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c return end c% function pggett (vpleft,vplevel,vphite,label,string) real vpleft ! left edge of window real vphite ! height of window real vplevel ! level of window character*(*) label ! optional label character*(*) string ! string to hold read c pggett ! returns right bound real xbox(4),ybox(4) character*1 cin integer pgcurs c if (vpleft.lt.0.0.or.vpleft.gt.1.0.or.vphite.lt.0.0.or.vphite.gt. $ 1.0) then write (*,*) 'PGGRABT: bad viewport value' return end if nchar=len(string) nlab=lenj(label) c call pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c chargoal=0.55 ! in world coords xtoffset=0.05 ! in world coords ytoffset=0.3 ! in world coords lens=lenj(string) if (lens.eq.0) string=' ' c vpright=vpleft+chargoal*float(nchar+nlab+1)*vphite vpright=amin1(vpright,1.0) c overide vpright=0.45 c overide pggett=vpright c call pgsvp (vpleft,vpright,vplevel-vphite/2.0,vplevel+vphite/2.0) call pgswin (0.0,1.0,0.0,1.0) call pgscf (2) call pgsls (1) call pgslw (2) call pgsch (1.0) call pgqtxt (0.0,0.0,0.0,0.0,'D',xbox,ybox) boxbase=amin1(ybox(1),ybox(2),ybox(3),ybox(4)) boxtop =amax1(ybox(1),ybox(2),ybox(3),ybox(4)) charhite=boxtop-boxbase charfact=chargoal/charhite call pgsch (charfact) c---loop getting values 30 icaught=0 10 call pgsci (1) call pgsfs (1) call pgrect (0.0,1.0,0.0,1.0) call pgsci (0) call pgstbg (1) if (nlab.gt.0) then call pgptxt (xtoffset,ytoffset,0.0,0.0,label) call pgqtxt (xtoffset,ytoffset,0.0,0.0,label,xbox,ybox) xlab=amax1(xbox(1),xbox(2),xbox(3),xbox(4)) call pgmove (xlab+xtoffset/2.0,0.0) call pgslw (9) call pgdraw (xlab+xtoffset/2.0,1.0) call pgslw (2) call pgqpos (xlab,ylab) xgof=xlab+xtoffset/2.0 else xgof=xtoffset end if xcurs=xgof ycurs=0.5 call pgptxt (xgof,ytoffset,0.0,0.0,string) ido=pgcurs (xcurs,ycurs,cin) if (xcurs.gt.0.0.and.xcurs.lt.(xgof-xtoffset/2.0).and.ycurs.lt.1.0 $ .and.ycurs.gt.0.0) then string=' ' go to 30 end if ihave=ichar(cin) if (ihave.eq.8.or.ihave.eq.127) then if (icaught.eq.0) icaught=1 string(icaught:nchar)=' ' icaught=icaught-1 if (icaught.lt.0) icaught=0 else if (ihave.eq.13) then go to 20 else icaught=icaught+1 string(icaught:icaught)=cin end if go to 10 c---loop getting values c 20 call pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c return end c% function pgputt (vpleft,vplevel,vphite,label,string) real vpleft ! left edge of window real vphite ! height of window real vplevel ! level of window character*(*) label ! optional label character*(*) string ! string holding output c pgputt ! returns right bound real xbox(4),ybox(4) c if (vpleft.lt.0.0.or.vpleft.gt.1.0.or.vphite.lt.0.0.or.vphite.gt. $ 1.0) then write (*,*) 'PGPUTT: bad viewport value' return end if nchar=len(string) nlab=lenj(label) c call pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c chargoal=0.55 ! in world coords xtoffset=0.05 ! in world coords ytoffset=0.3 ! in world coords c vpright=vpleft+chargoal*float(nchar+nlab+1)*vphite vpright=amin1(vpright,1.0) c override vpright=0.95 c override pgputt=vpright c call pgsvp (vpleft,vpright,vplevel-vphite/2.0,vplevel+vphite/2.0) call pgswin (0.0,1.0,0.0,1.0) call pgscf (2) call pgsls (1) call pgslw (2) call pgsch (1.0) call pgqtxt (0.0,0.0,0.0,0.0,'D',xbox,ybox) boxbase=amin1(ybox(1),ybox(2),ybox(3),ybox(4)) boxtop =amax1(ybox(1),ybox(2),ybox(3),ybox(4)) charhite=boxtop-boxbase charfact=chargoal/charhite call pgsch (charfact) call pgsci (1) call pgsfs (1) call pgrect (0.0,1.0,0.0,1.0) call pgsci (0) call pgstbg (1) if (nlab.gt.0) then call pgptxt (xtoffset,ytoffset,0.0,0.0,label) call pgqtxt (xtoffset,ytoffset,0.0,0.0,label,xbox,ybox) xlab=amax1(xbox(1),xbox(2),xbox(3),xbox(4)) call pgmove (xlab+xtoffset/2.0,0.0) call pgslw (9) call pgdraw (xlab+xtoffset/2.0,1.0) call pgslw (2) call pgqpos (xlab,ylab) xgof=xlab+xtoffset/2.0 else xgof=xtoffset end if call pgptxt (xgof,ytoffset,0.0,0.0,string) c call pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c return end c% subroutine pgbutx (padllx,padlly,padurx,padury) real padllx,padlly,padurx,padury ! BB of button pad c if (padllx.lt.0.0.or.padurx.gt.1.0.or.padlly.lt.0.0.or.padury.gt. $ 1.0) then write (*,*) 'PGBUTX: bad viewport values' return end if c call pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c call pgsvp (padllx,padurx,padlly,padury) call pgswin (0.0,1.0,0.0,1.0) call pgsci (0) call pgsfs (1) call pgrect (0.0,1.0,0.0,1.0) c c call pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c return end c% subroutine pgdiax (vpleft,vplevel,vphite,vpright) real vpleft ! left edge of window real vphite ! height of window real vplevel ! level of window real vpright ! right edge from pggett etc. c if (vpleft.lt.0.0.or.vpleft.gt.1.0.or.vphite.lt.0.0.or.vphite.gt. $ 1.0.or.vpright.lt.0.0.or.vpright.gt.1.0) then write (*,*) 'PGDIAX: bad viewport values' return end if c call pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c call pgsvp (vpleft,vpright,vplevel-vphite/2.0,vplevel+vphite/2.0) call pgswin (0.0,1.0,0.0,1.0) call pgsci (0) call pgsfs (1) call pgrect (0.0,1.0,0.0,1.0) c c call pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) c return end c% subroutine pgstord (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) call pgqcf (igfont) call pgqch (ghite) call pgqci (igcolr) call pgqtbg (igtxtc) call pgqfs (igfils) call pgqls (iglins) call pgqlw (iglinw) call pgqvp (0,vpllx,vpurx,vplly,vpury) call pgqwin (wllx,wurx,wlly,wury) !needs to be set? call pgqpos (xpos,ypos) !needs to be set? return end c% subroutine pgputbd (igfont,ghite,igcolr,igtxtc,igfils,iglins, $ iglinw,vpllx,vpurx,vplly,vpury,wllx,wurx,wlly,wury,xpos,ypos) call pgscf (igfont) call pgsch (ghite) call pgsci (igcolr) call pgstbg (igtxtc) call pgsfs (igfils) call pgsls (iglins) call pgslw (iglinw) call pgsvp (vpllx,vpurx,vplly,vpury) call pgswin (wllx,wurx,wlly,wury) call pgmove (xpos,ypos) return end c% integer function lenj(string) character*(*) string c max=len(string) do i=max,1,-1 if (ichar(string(i:i)).ge.33.and.ichar(string(i:i)).le.126) x then lenj=i return end if end do lenj=0 return end