5     ! ####################################################
10    !           Main
15    ! Change normalization scheme by DATA in Rmnsb, Rmnjf & Presphe.
20    !  ?/2106
25    ! ####################################################
30    Symbol(0)
35    PRINTER IS CRT
40    OUTPUT KBD;CHR$(255)&CHR$(75);   ! Clear screen.
45    Ct=3
50    Cl=5
55    Cr=45
60  Menu: PRINT TABXY(20,1),"*** Program Menu ***"
65    ! ------------------Scattering--------------------------
70    L=2
75    PRINT TABXY(Ct,L),"Scattering"
80    PRINT TABXY(Cl,L+1),"LP:L Pattern (Scatpatl)"
85    PRINT TABXY(Cl,L+2),"P3:SLR Pattern (Scatpat)"
90    PRINT TABXY(Cl,L+3),"K3:SLR-ka (Scatka)"
95    PRINT TABXY(Cl,L+4),"L1:L 1dir (Scatl1d)"
100   PRINT TABXY(Cl,L+5),"DF:L Dens (Densf)"
105   PRINT TABXY(Cl,L+6),"BB:Body,Blad.-kL/2 (Bodyblad)"
110   PRINT TABXY(Cl,L+7),"SP:Seq.Back.Pat.(Seqbackpat)"
115   PRINT TABXY(Cr,L+1),"TM:NOF ave. meas. TS (Tsavemeas)"
120   PRINT TABXY(Cr,L+2),"TF:A-L/Lamda (Tsavef)"
125   ! ------------------SWF---------------------------------
130   L=10
135   PRINT TABXY(Ct,L),"SWF"
140   PRINT TABXY(Cl,L+1),"RH:Rmnhc"
145   PRINT TABXY(Cl,L+2),"DC:D_check"
150   PRINT TABXY(Cl,L+3),"RC:Rmnc"
155   PRINT TABXY(Cl,L+4),"SR: Smn & Rmn - ka"
160   PRINT TABXY(Cl,L+5),"LC: Lamda check "
165   PRINT TABXY(Cl,L+6),"LF: Make or Store Lamda file"
170   ! ------------------Other-------------------------------
175   L=17
180   PRINT TABXY(Ct,L),"Other"
185   PRINT TABXY(Cl,L+1),"PF:Plot by filed data"
190   ! ------------------Select pro.-------------------------
195 ! INPUT "Which FN for Lamda (F/C)?",Lcf$
200   DISP "Select pro. by cord. Cord name ";
205   LINPUT Pro$
210   OUTPUT KBD;CHR$(255)&CHR$(75);   ! Clear screen.
215  More: SELECT Pro$
220   CASE "LP"
225     Scatpatl
230   CASE "P3"
235     Scatpat
240   CASE "K3"
245     Scatka
250   CASE "L1"
255     Scatl1d
260   CASE "DF"
265     Densf
270   CASE "BB"
275     Bodyblad
280   CASE "RH"
285     Rmnhc
290   CASE "RC"
295     Rmn_c2
300   CASE "DC"
305     D_c
310   CASE "SR"
315     Srmnc
320   CASE "PF"
325     Pltfile(0)
330   CASE "LC"
335     Lamdac
340   CASE "LF"
345     Lfile
350   CASE "SP"
355     Seqbackpat
360   CASE "TM"
365     Tsavemeas
370   CASE "TF"
375     Tsavef
380   CASE ELSE
385     GOTO Menu
390   END SELECT
395   DISP "Another pro. to <CONT>"
400   PAUSE
405   GOTO 40
410   END
415   ! ####################################################
420   DEF FNSmn(M,N,C,Eta)
425   ! Angle function Smn.  See Flammer eq.(3.1.3a).
430   ! #####################(S)############################
435     COM /D/Dd(-1000:1000),Rmin,Rmid,Rstart,Rmax   ! <--Dbyd2,Dnorm,Dminus
440     IF Eta=0 THEN Eta=1.E-5
445     IF C=0 THEN RETURN FNLgna1(M,N,0,Eta)
450     S=0
455     FOR R=Rstart TO Rmax STEP 2
460       S=S+Dd(R)*FNLgna1(M,M+R,0,Eta)
465     NEXT R
470     RETURN S
475   FNEND
480   ! ####################################################
485   SUB Presphe(M,N,C,L)
490     !  Calculate Lamda & d, and gives Rstart & Rmax in COM below.
495     !  Lpre + Dbyd2 + Dnorm.
500     !  L must be return for calculation of d(-r).
505     !  FNLamda or FNLamdaf is selected by COM /Lcf/.
510     ! ###################(T)##############################
515     COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
520     COM /Lcf/Lcf$[8]
525     READ Norm$
530     DATA "F"
535     IF Lcf$="C" THEN CALL Lpre(M,N)
540     Dbyd2(M,N,C,L)
545     IF Norm$="F" THEN
550       Dnorm(M,N)
555     ELSE
560       Dnorm2(M,N)
565     END IF
570   SUBEND
575   ! ####################################################
580   DEF FNSnorm(M,N)
585   ! Norm of Smn. See Flammer eq.(3.1.33).
590   ! #####################(S2)###########################
595     COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
600     M2=2*M
605     Nume=FNFact(Rstart+M2)
610     Deno=Rstart*2+M2+1
615     Def=1
620     S=Nume/Deno*Dd(Rstart)^2
625     FOR R=Rstart+2 TO Rmax STEP 2
630       Nume=Nume*(R+M2)*(R-1+M2)
635       Deno=Deno+4
640       Def=Def*R*(R-1)
645       S=S+Nume/Deno/Def*Dd(R)^2
650     NEXT R
655     RETURN S*2
660   FNEND
665   ! ####################################################
670   SUB Rmnh(Kind,M,N,C,L,Xi,Rmn1,Rmn2,Rmn1d,Rmn2d,Err)
675   ! Spheroidal radial functions combined according to Hanish et al.
680   ! See Hanish et al. p.xvi.
685   ! Rmn1 always by spherical Bessel function.
690   ! Rmn2 are selected by C & Xi.
695   ! Presphe must be called beforehand.
700   ! Err is relative Wronskian error in %.  =0 for Kind=1,10.
705   ! ------------------------------------------------------
710   ! Kind   Rmn1  Rmn2  Rmn1d  Rmn2d
715   !   1      O
720   !   2      O     O      O      O
725   !   3      O     O      O      O
730   !   10     O            O
735   !   20     O     O      O      O
740   !   30     O     O      O      O
745   ! ####################################################
750     DISP "Now in Rmnh."
755     COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
760     IF Kind=1 OR Kind=10 THEN
765       Rmnsb(M,N,C,Xi,Kind,Rmn1,D,Rmn1d,D)
770       Err=0
775       SUBEXIT
780     END IF
785     SELECT Xi
790     CASE >5
795       IF C>.1 THEN
800         GOSUB Sb
805       ELSE
810         GOSUB Jf
815       END IF
820     CASE >2
825       IF C>1 THEN
830         GOSUB Sb
835       ELSE
840         GOSUB Jf
845       END IF
850     CASE >1.3
855       IF C>10 THEN
860         GOSUB Sb
865       ELSE
870         GOSUB Jf
875       END IF
880     CASE >1.05
885       GOSUB Sbjf
890     CASE >=1
895       GOSUB Jf
900     CASE ELSE
905       BEEP
910       PRINT "Input error in Rmnh."
915     END SELECT
920     IF ABS(Err)>5 THEN
925       BEEP
930       BEEP
935       PRINT "Wronskian error ";Err;"% in Rmnh. m=";M;" n=";N;" h=";C;" Xi=";Xi
940     END IF
945     SUBEXIT
950  Sb: ! ------------------------------------------------------
955     Rmnsb(M,N,C,Xi,30,Rmn1,Rmn2,Rmn1d,Rmn2d)
960     Err=FNWronerr(C,Xi,Rmn1,Rmn2,Rmn1d,Rmn2d)
965     Js$="S"                 !**************
970     RETURN
975  Jf: ! ------------------------------------------------------
980     Rmnsb(M,N,C,Xi,10,Rmn1,D,Rmn1d,D)    ! Rmn1 by Sb.
985     Dminus(M,N,C,L)
990     Rmnjf(M,N,C,Xi,20,D,D,Rmn2,D,Rmn2d)  ! Rmn2 by Jf.
995     Err=FNWronerr(C,Xi,Rmn1,Rmn2,Rmn1d,Rmn2d)
1000    Js$="J"                 !*****************
1005    RETURN
1010  Sbjf:! ------------------Compare both method-----------------
1015    Rmnsb(M,N,C,Xi,30,Rmn1,Rmn2s,Rmn1d,Rmn2ds)
1020    Errs=FNWronerr(C,Xi,Rmn1,Rmn2s,Rmn1d,Rmn2ds)
1025    Dminus(M,N,C,L)
1030    Rmnjf(M,N,C,Xi,20,D1,D2,Rmn2j,D3,Rmn2dj)
1035    Errj=FNWronerr(C,Xi,Rmn1,Rmn2j,Rmn1d,Rmn2dj)
1040  ! PRINT "Rmn2jf in Rmnh.",Rmn2j  !**********
1045    IF ABS(Errj)<ABS(Errs) THEN         ! Jf.
1050      Rmn2=Rmn2j
1055      Rmn2d=Rmn2dj
1060      Err=Errj
1065      Js$="J"              !***********
1070    ELSE                                ! Sb.
1075      Rmn2=Rmn2s
1080      Rmn2d=Rmn2ds
1085      Err=Errs
1090      Js$="S"              !*************
1095    END IF
1100    RETURN
1105  SUBEND
1110     ! ####################################################
1115  DEF FNWronerr(C,Xi,Rmn1,Rmn2,Rmn1d,Rmn2d)
1120     ! Relative Wronskian error in %.
1125     ! ####################################################
1130    Cal=Rmn1*Rmn2d-Rmn2*Rmn1d
1135    Theo=1/C/(Xi*Xi-1)
1140    RETURN (Cal-Theo)/Theo*100
1145  FNEND
1150     ! ####################################################
1155  SUB Lpre(M,N)
1160  !  Expansion coef. of lamda for small c.
1165  !  Preparation for FNLamdasc.  See Flammer p.18-19.
1170  !  This SUB must be called before FNLamdasc when m & n are changed.
1175  !  Need FNBc
1180  !  L for large 2k greatly affects Lamda for large c.
1185  ! #####################(L)############################
1190    COM /Lv/Lv(10)           ! Common with FNLamdasc.
1195    DIM Nn(-9:11),Nmm(-3:4),Nmp(-3:4)
1200    FOR I=-9 TO 11 STEP 2
1205      Nn(I)=2*N+I            ! Nn=2n-9 TO 2n+1
1210    NEXT I
1215    FOR I=-3 TO 4
1220      Nmm(I)=N-M+I           ! Nmm=n-m-3 TO n-m+4
1225    NEXT I
1230    FOR I=-3 TO 4
1235      Nmp(I)=N+M+I           ! Nmp=n+m-3 TO n+m+4
1240    NEXT I
1245    M4=4*M^2-1               ! M4=4m^2-1
1250    ! ------------------------------------------------------
1255    Lv(0)=N*(N+1)
1260    Lv(2)=(1-M4/Nn(-1)/Nn(3))/2
1265    Lv(4)=(Nmm(-1)*Nmm(0)*Nmp(-1)*Nmp(0)/Nn(-3)/Nn(-1)^3/Nn(1)-Nmm(1)*Nmm(2)*Nmp(1)*Nmp(2)/Nn(1)/Nn(3)^3/Nn(5))/2
1270    Lv(6)=M4*(Nmm(1)*Nmm(2)*Nmp(1)*Nmp(2)/Nn(-1)/Nn(1)/Nn(3)^5/Nn(5)/Nn(7)-Nmm(-1)*Nmm(0)*Nmp(-1)*Nmp(0)/Nn(-5)/Nn(-3)/Nn(-1)^5/Nn(1)/Nn(3))
1275    A=Nmm(-1)*Nmm(0)*Nmp(-1)*Nmp(0)/Nn(-5)^2/Nn(-3)/Nn(-1)^7/Nn(1)/Nn(3)^2-Nmm(1)*Nmm(2)*Nmp(1)*Nmp(2)/Nn(-1)^2/Nn(1)/Nn(3)^7/Nn(5)/Nn(7)^2
1280    B=Nmm(-3)*Nmm(-2)*Nmm(-1)*Nmm(0)*Nmp(-3)*Nmp(-2)*Nmp(-1)*Nmp(0)/Nn(-7)/Nn(-5)^2/Nn(-3)^3/Nn(-1)^4/Nn(1)
1285    C=Nmm(1)*Nmm(2)*Nmm(3)*Nmm(4)*Nmp(1)*Nmp(2)*Nmp(3)*Nmp(4)/Nn(1)/Nn(3)^4/Nn(5)^3/Nn(7)^2/Nn(9)
1290    D=Nmm(1)^2*Nmm(2)^2*Nmp(1)^2*Nmp(2)^2/Nn(1)^2/Nn(3)^7/Nn(5)^2-Nmm(-1)^2*Nmm(0)^2*Nmp(-1)^2*Nmp(0)^2/Nn(-3)^2/Nn(-1)^7/Nn(1)^2
1295    E=Nmm(-1)*Nmm(0)*Nmm(1)*Nmm(2)*Nmp(-1)*Nmp(0)*Nmp(1)*Nmp(2)/Nn(-3)/Nn(-1)^4/Nn(1)^2/Nn(3)^4/Nn(5)
1300    Lv(8)=2*M4^2*A+(B-C)/16+D/8+E/2
1305    A=Nn(-5)*Nn(-1)^2*Nn(3)/M4*Lv(6)-4*Lv(4)+(6*N-19)*FNBc(M,N-M-2)/2/Nn(-3)/Nn(-9)+16*M4^2/Nn(-5)^2/Nn(-1)^3/Nn(3)^2     ! Only case of R<2.
1310    B=Nn(-1)*Nn(3)^2*Nn(7)/M4*Lv(6)-4*Lv(4)-(6*N-25)*FNBc(M,N-M+4)/2/Nn(5)/Nn(11)-16*M4^2/Nn(-1)^2/Nn(3)^3/Nn(7)^2
1315    Lv(10)=-M4*FNBc(M,N-M)/4/Nn(-5)/Nn(-1)^4/Nn(3)*A-M4*FNBc(M,N-M+2)/4/Nn(-1)/Nn(3)^4/Nn(7)*B
1320         ! L10 is not fully expanded.
1325  SUBEND
1330  ! ####################################################
1335  DEF FNBc(M,R)
1340  !  Beta/c^4. See Flammer eq.(3.1.6).
1345  !  There is CSUB alternative Pbc.
1350  ! #####################(C)############################
1355    M2r=2*M+R
1360    M2r2=2*M+2*R
1365    RETURN R*(R-1)*M2r*(M2r-1)/(M2r2-1)^2/(M2r2-3)/(M2r2+1)
1370  FNEND
1375  ! ####################################################
1380  SUB Dnorm(M,N)
1385  !  Normalization of d/d.  See Flammer 3.1.5.
1390  ! #####################(D)############################
1395    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
1400    DISP "Now in Dnorm."
1405    SELECT Rstart
1410    CASE 0
1415      A=1
1420      IF M=0 THEN GOTO Skip1
1425      FOR I=M+1 TO 2*M
1430        A=A*I
1435      NEXT I
1440    CASE 1
1445      A=1/2
1450      FOR I=M+2 TO 2*M+2
1455        A=A*I
1460      NEXT I
1465    END SELECT
1470  Skip1: IF Rstart=N-M THEN B=A
1475    Sadd=A*Dd(Rstart)
1480    FOR R=Rstart+2 TO Rmax STEP 2
1485      SELECT Rstart
1490      CASE 0
1495        A=-A*(2*M+R-1)*(2*M+R)/2/R/(M+R/2)
1500      CASE 1
1505        A=-A*(2*M+R+1)*(2*M+R)/2/(R-1)/(M+(R+1)/2)
1510      END SELECT
1515      IF R=N-M THEN B=A
1520      Sadd=Sadd+A*Dd(R)
1525    NEXT R
1530    Dmn=B/Sadd
1535    FOR R=Rstart TO Rmax STEP 2
1540      Dd(R)=Dmn*Dd(R)
1545    NEXT R
1550  SUBEND
1555  ! ####################################################
1560  SUB Rmnsb(M,N,C,Xi,K,Rmn1,Rmn2,Rmn1d,Rmn2d)
1565  !  Spheroidal radial wave functions & their derivatives                   6415  !    by spherical Bessel functions.
1570  !  See Flammer eqs.(4.1.15) & (4.1.19).
1575  !  K=1      Rmn1
1580  !    2      Rmn2
1585  !    3      Rmn1, Rmn2
1590  !    10     Rmn1, Rmn1d
1595  !    20     Rmn2, Rmn2d
1600  !    30     Rmn1, Rmn2, Rmn1d, Rmn2d
1605  !  Dd,Rmax,Rstart must be given by COM.
1610  !  Needs  SUBSbes_der, FNSjbes, Snbes, & Fact.
1615  ! #####################(R)############################
1620    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
1625    READ Norm$
1630    DATA "F"
1635    DISP "Now in Rmnsb."
1640  ! ------------------Initial value-----------------------
1645    SELECT Rstart
1650    CASE 0                        ! N-M:even
1655      Sign=(-1)^((N-M)/2)         ! Initial value is for Rstart.
1660      Mm=FNFact(2*M)
1665    CASE 1                        ! N-M:odd
1670      Sign=(-1)^((N-M-1)/2)
1675      Mm=FNFact(2*M+1)
1680    END SELECT
1685    Cxi=C*Xi
1690    Xx=1-1/Xi/Xi
1695    Rmn1=0
1700    Rmn2=0
1705    Rmn1d=0
1710    Rmn2d=0
1715    S1=0                          ! For 1st kind.
1720    S2=0                          ! For 2nd kid.
1725    S1d=0                         ! For der. of 1st.
1730    S2d=0                         ! For der. of 2nd.
1735    Sdmr=0                        ! For denominator.
1740  ! ------------------Sumation----------------------------
1745    FOR R=Rstart TO Rmax STEP 2
1750      SELECT R
1755      CASE Rstart
1760        Rr=1
1765      CASE ELSE
1770        Rr=Rr*R*(R-1)
1775        Sign=Sign*(-1)
1780        Mr=2*M+R
1785        Mm=Mm*Mr*(Mr-1)
1790      END SELECT
1795      Dmr=Dd(R)*Mm/Rr
1800      Sdmr=Sdmr+Dmr
1805  !  PRINT "M,N,R,Sdmr,Rmax in Rmnsb.",M,N,R,PROUND(Sdmr,-2),Rmax  !************
1810      Sind=Sign*Dmr
1815      IF K=1 OR K=3 THEN S1=S1+Sind*FNSjbes(M+R,Cxi)
1820      IF K=2 OR K=3 THEN S2=S2+Sind*FNSnbes(M+R,Cxi)
1825      IF K=10 OR K=30 THEN
1830        Sbes_der(M+R,Cxi,"J",1,Sj,Sjd,Sjdd)
1835        S1=S1+Sind*Sj
1840        S1d=S1d+Sind*C*Sjd
1845      END IF
1850      IF K=20 OR K=30 THEN
1855        Sbes_der(M+R,Cxi,"N",1,Sn,Snd,Sndd)
1860        S2=S2+Sind*Sn
1865        S2d=S2d+Sind*C*Snd
1870      END IF
1875    NEXT R
1880  ! ------------------Result------------------------------
1885    IF Norm$="F" THEN
1890      A=Xx^(M/2)/Sdmr
1895    ELSE
1900      A=Xx^(M/2)*FNFact(N-M)/FNFact(N+M)!************
1905    END IF
1910    IF K>=10 THEN B=A/Xx*M/Xi^3
1915    IF K<>2 AND K<>20 THEN Rmn1=A*S1
1920    IF K<>1 AND K<>10 THEN Rmn2=A*S2
1925    IF K=10 OR K=30 THEN Rmn1d=A*S1d+B*S1
1930    IF K=20 OR K=30 THEN Rmn2d=A*S2d+B*S2
1935  SUBEND
1940  ! ####################################################
1945  SUB Dminus(M,N,C,L)
1950  !  d for -r & d(rho|r).  See Flammer 3.5.
1955  !  Rlast < Rmin < Rmid < Rstart < Rmax < Rend
1960  !           |d(rho,r)| d(-r)| d(r) |
1965  ! #####################(D)############################
1970    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
1975    DISP "Now in Dminus."
1980    DATA 1E-10,1E-2,1E-10
1985    READ Ddmin,Cuterror,Rho  ! Small Cuterror waste time.
1990  ! INPUT "Rho",Rho    !***********
1995    ! ------------------d(-r)/d(-r+2)-----------------------
2000    K=-(2*M-Rstart)          ! Last term of d(-r)
2005    Rmid=K
2010    IF M=0 THEN GOTO Rho     ! No d(-r).
2015    Dd(K)=FNDdlast(M,L,K,C)  ! Dbyd for Rmid.
2020    IF M=1 THEN GOTO D       ! Only one d(-r).
2025    REPEAT
2030      K=K+2                  ! K increases from Rmid+2 to Rstart-2.
2035      Kp2=K+2
2040      Dd(K)=FNDdr(M,L,K,C,Dd(K-2))  ! d(-r)/d(-r+2)
2045    UNTIL Kp2=Rstart
2050    ! ------------------d/d to d----------------------------
2055  D: FOR R=Rstart-2 TO Rmid STEP -2
2060      Dd(R)=Dd(R)*Dd(R+2)
2065    NEXT R
2070    ! ------------------Determin Rmin-----------------------
2075  Rho: Rmin=Rmid-6              ! Provisional.
2080    Rlast=Rmin-4                ! Do.
2085  Again: GOSUB Ddr
2090    GOSUB D2
2095    IF ABS(Dd(Rmin))>Ddmin THEN
2100      Rmin=Rmin-4
2105      Rlast=Rmin-4
2110      DISP "Rmin in Dminus.",Rmin  !**********
2115      GOTO Again
2120    END IF
2125    Rlast0=Rlast
2130    ! ------------------Determin Rlast----------------------
2135    REPEAT                 ! Until 2nd term of deno. of eq.(3.5.6) is
2140      Rlast=Rlast-2        !   largely less than 1st term.
2145      DISP "Rlast in Dminus.",Rlast  !********
2150    UNTIL ABS(FNLc(M,Rlast,C)*FNDdlast(M,L,Rlast,C)/FNLb(M,L,Rlast+2,C))<Cuterror
2155    ! ------------------Calculate d(rho|r)------------------
2160    IF Rlast=Rlast0 THEN SUBEXIT    ! Already calculated.
2165    GOSUB Ddr
2170    GOSUB D2
2175    SUBEXIT
2180    ! ----------------d(-r)/d(-r+2)-------------------------
2185  Ddr: Dd(Rlast)=FNDdlast(M,L,Rlast,C)! Cut off.
2190    FOR R=Rlast+2 TO Rmid-4 STEP 2
2195      Dd(R)=FNDdr(M,L,R,C,Dd(R-2))
2200    NEXT R
2205    Dd(Rmid-2)=-FNNd(M,Rmid+Rho,C)/(FNLb(M,L,Rmid-2,C)+FNLc(M,Rmid-4,C)*Dd(Rmid-4))/Rho
2210    RETURN
2215    ! ----------------d/d to d(rho/r)-----------------------
2220  D2: FOR R=Rmid-2 TO Rmin STEP -2
2225      Dd(R)=Dd(R)*Dd(R+2)
2230    NEXT R
2235    RETURN
2240  SUBEND
2245  ! ####################################################
2250  DEF FNDdlast(M,L,R,C)
2255  !  d/d for Rmid or Rlast.  See Flammer eq.(3.5.6).
2260  !  -A(-2m+2(+3))/B(-2m(+1)).
2265  ! #####################(D)############################
2270    RETURN -FNNd(M,R+2,C)/FNLb(M,L,R,C)
2275  FNEND
2280  ! ####################################################
2285  DEF FNDdr(M,L,R,C,Ddprev)
2290  !  d/d for -r.  d(-r)/d(-r+2).  See Flammer eq.(3.5.6).
2295  !  Ddprev is previous d/d.  Nd equales A.
2300  ! #####################(D)############################
2305    RETURN -FNNd(M,R+2,C)/(FNLb(M,L,R,C)+FNLc(M,R-2,C)*Ddprev)
2310  FNEND
2315  ! ####################################################
2320  DEF FNLb(M,L,R,C)
2325  !  Coef. of difference eq. of d, B.  See Flammer eq.(3.5.4).
2330  ! #####################(C)############################
2335    Mr=M+R
2340    Mr2=2*Mr
2345    RETURN Mr*(Mr+1)-L+(Mr2*(Mr+1)-2*M*M-1)/(Mr2-1)/(Mr2+3)*C*C
2350  FNEND
2355  ! ####################################################
2360  DEF FNLc(M,R,C)
2365  !  Coef. of difference eq. of d, C.  See Flammer eq.(3.5.5).
2370  ! #####################(C)############################
2375    Mr=2*(M+R)
2380    RETURN (R+1)*(R+2)/(Mr+1)/(Mr+3)*C*C
2385  FNEND
2390  ! ####################################################
2395  SUB Rmnjf(M,N,C,Xi,K,Smn,Rmn1,Rmn2,Rmn1d,Rmn2d)
2400  !  Spheroidal radial wave functions & derivatives by joining factors.
2405  !  See Flammer eqs.(4.2.1),(4.2.4),& (4.2.6).
2410  !  Smn is angle function which is used to get Rmn1 (Result).              8925  !  K is same as in SUBRmnsb.
2415  !  SUBDminus must be called beforehand to get Rmn2 & Rmn2d.
2420  ! #####################(R)############################
2425    DISP "In Rmnjf"
2430    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
2435    DIM Q(20,-20:20),Qd(20,-20:20)
2440    READ Norm$
2445    DATA "F"
2450    Smn=0
2455    Rmn1=0
2460    Rmn2=0
2465    Rmn1d=0
2470    Rmn2d=0
2475    ! ------------------Kappa-------------------------------
2480    IF Norm$="F" THEN
2485      Kappa(M,N,C,3,Kap1,Kap2)
2490    ELSE
2495      Kappa(M,N,C,1,Kap1,Kap2)
2500      Kappa2(M,N,C,Kap2)            ! Kappa2 by Hanish et al.
2505    END IF
2510    ! ------------------Rmn1--------------------------------
2515    IF K=1 OR K=3 THEN      ! First kind.
2520      Smn=FNSmn(M,N,C,Xi)
2525      Rmn1=Smn/Kap1
2530    END IF
2535    ! ------------------Rmn2--------------------------------
2540    IF K=2 OR K=3 THEN      ! Second kind.
2545      Lgna2(M,M+Rmax,(-M+Rstart)*(M<>0),Xi,Q(*))
2550      !               When M=0, negative n is not needed.
2555      S1=0
2560      FOR R=Rmid TO Rmax STEP 2
2565        S1=S1+Dd(R)*Q(M,M+R)     ! Need Q(M,-N)
2570      NEXT R
2575      S2=0
2580      FOR R=Rmid-2 TO Rmin STEP -2
2585        S2=S2+Dd(R)*FNLgna1(M,-R-M-1,0,Xi)
2590      NEXT R
2595  !   PRINT "Rmid,Kap2,Q-1,Q1,S1,S2",Rmid,Kap2,Q(1,-1),Q(1,1),S1,S2  !******
2600      Rmn2=(S1+S2)/Kap2
2605    END IF
2610    ! ------------------Rmn1 & Rmn1'------------------------
2615    IF K=10 OR K=30 THEN          ! 1st kind & derivative.
2620      Smnd(M,N,C,Xi,Smn,Smnd)
2625      Rmn1=Smn/Kap1
2630      Rmn1d=Smnd/Kap1
2635    END IF
2640    ! ------------------Rmn2 & Rmn2'------------------------
2645    IF K=20 OR K=30 THEN          ! 2nd kind & derivative.
2650      S1=0
2655      S1d=0
2660      Lgna2_der(M,M+Rmax,MIN((M+Rmid),0),Xi,Q(*),Qd(*))
2665      FOR R=Rmid TO Rmax STEP 2
2670   !    CALL Lgna_der(M,M+R,Xi,"Q",Qa,Qd)
2675        S1=S1+Dd(R)*Q(M,M+R)
2680        S1d=S1d+Dd(R)*Qd(M,M+R)
2685   !    PRINT USING """R,d,Q,Qd"",2X,S2D,2X,3(SD.4DE,2X)";R,Dd(R),Qa,Dd(R)*Qa                 !**************
2690      NEXT R
2695      S2=0
2700      S2d=0
2705      FOR R=Rmid-2 TO Rmin STEP -2
2710        Lgna_der(M,-R-M-1,Xi,"P",Pa,Pd)
2715        S2=S2+Dd(R)*Pa
2720        S2d=S2d+Dd(R)*Pd
2725   !    PRINT USING """R,d,P,Pd"",2X,S2D,2X,3(SD.4DE,2X)";R,Dd(R),Pa,Dd(R)*Pa      !****************
2730      NEXT R
2735      Rmn2=(S1+S2)/Kap2
2740      Rmn2d=(S1d+S2d)/Kap2
2745    END IF
2750  SUBEND
2755  ! ####################################################
2760  SUB Smnd(M,N,C,Eta,Smn,Smnd)
2765  ! Smn & its derivative.
2770  ! #####################(S3)###########################
2775    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
2780    IF Eta=0 THEN Eta=1.E-10
2785    IF C=0 THEN
2790      Lgna_der(M,N,Eta,"P",P,Pd)
2795      Smn=P
2800      Smnd=Pd
2805      SUBEXIT
2810    END IF
2815    S=0
2820    Sd=0
2825    FOR R=Rstart TO Rmax STEP 2
2830      Lgna_der(M,M+R,Eta,"P",P,Pd)
2835      S=S+Dd(R)*P
2840      Sd=Sd+Dd(R)*Pd
2845    NEXT R
2850    Smn=S
2855    Smnd=Sd
2860  SUBEND
2865  ! ####################################################
2870  DEF FNLamda(M,N,C,Accuracy)
2875  !  Lamda(m,n,c). See Flammer eq.(3.1.17) & 3.1.4.
2880  !  Need FNL,Bc,Nmru,Nmrl,Lamdasc,Lamdalc & SUBRendfind.
2885  !  Accuracy is only transfered to Bouwkamp or Pincer.
2890  !         |-Lamdasc                                   C<CtlO
2895  !         |-Lamdasc--Pincer or Bouwkamp               Ctl0<C<Ctl
2900  !  Lamda----Intmed--Initial Lamda--Bouwkamp or Pincer Ctl<C<Cth
2905  !         |-Lamdalc--Pincer or Bouwkamp               Cth<C<Cth0
2910  !         |-Lamdalc                                   C>Cth0
2915  ! #####################(L)############################
2920    COM /C/C4
2925    DATA .5,2,20,100,10,B
2930    READ Ctl0,Ctl,Cth,Cth0,Stepn,Bp$  ! Bp$:B for Bouwkamp, P for Pincer.
2935    C4=C^4
2940    IF M<0 OR N<M THEN
2945      BEEP
2950      OUTPUT 1;"m>=0 and n>=m must be satisfied in Lamda. Return 0."
2955      RETURN 0
2960    END IF
2965    ! ------------------------------------------------------
2970    SELECT C
2975    CASE <Ctl0                    ! Only by power series expansion of c.
2980      RETURN FNLamdasc(M,N,C)
2985    CASE Ctl0 TO Ctl
2990      IF Bp$="P" THEN
2995        RETURN FNPincer(M,N,C,FNLamdasc(M,N,C),Accuracy)
3000      ELSE
3005        RETURN FNBouwkamp(M,N,C,FNLamdasc(M,N,C),Accuracy)
3010      END IF
3015    CASE Ctl TO Cth               ! Select Ll or Lh & Bowukamp's or pincer
3020                                  !   method.
3025      GOSUB Intmed
3030      RETURN L
3035    CASE Cth TO Cth0
3040      IF Bp$="P" THEN
3045        RETURN FNPincer(M,N,C,FNLamdalc(M,N,C),Accuracy)
3050      ELSE
3055        RETURN FNBouwkamp(M,N,C,FNLamdalc(M,N,C),Accuracy)
3060      END IF
3065    CASE >Cth0                    ! Only by power series expansion of 1/c.
3070      RETURN FNLamdalc(M,N,C)
3075    END SELECT
3080    ! ------------------------------------------------------
3085  Intmed:  ! Find C which gives minimum of ABS(Lamdasc-Lamdalc),Cmin.
3090    DISP "Now in FNLamda,Intmed."
3095    Cinc=(Cth/Ctl)^(1/Stepn)
3100    Cc=Ctl
3105    Dcmin=10000
3110    WHILE Cc<=Cth
3115      Dc=ABS(FNLamdasc(M,N,Cc)-FNLamdalc(M,N,Cc))
3120      IF Dc<Dcmin THEN
3125        Dcmin=Dc
3130        Cmin=Cc                   ! Result.
3135      END IF
3140      Cc=Cc*Cinc
3145    END WHILE
3150    SELECT C                      ! Initial Lamda.
3155    CASE <Cmin
3160      Linit=FNLamdasc(M,N,C)
3165    CASE >=Cmin
3170      Linit=FNLamdalc(M,N,C)
3175    END SELECT
3180    IF Bp$="B" THEN L=FNBouwkamp(M,N,C,Linit,Accuracy)
3185    IF Bp$="P" THEN L=FNPincer(M,N,C,Linit,Accuracy)
3190    RETURN
3195  FNEND
3200  ! ####################################################
3205  DEF FNNd(M,R,C)
3210  !  Coef. of difference eq. of d.
3215  !  N/(d/d) or A. See Flammer eq.(3.1.7) or (3.5.3).
3220  ! #####################(D)############################
3225    Mr=2*M+R
3230    Mr2=Mr+R
3235    RETURN Mr*(Mr-1)/(Mr2-1)/(Mr2+1)*C^2
3240  FNEND
3245  ! ####################################################
3250  DEF FNPincer(M,N,C,Linit,Accu)
3255  ! Pincer method for Lamda.
3260  ! Return Lamda.
3265  ! #####################(L)############################
3270    DIM U(-1:1),L(-1:1)
3275    DATA 5E-4,10    ! Large & small Initaccu waste time.
3280    READ Initaccu,Accufac
3285      ! Initaccu must be sufficiently small. U has some roots.
3290    Times=0                       ! Times of iteration of Tloop.
3295    Rendp=0                       ! Previous Rend.
3300    Lmd=Linit
3305  ! ------------------------------------------------------
3310  Tloop: Rendfind(M,N-M+2,C,Lmd,Rend)
3315    IF Rend=Rendp THEN            ! Exit here.
3320 !  OUTPUT 1 USING Fmt;Lmd,Times,Count,Curaccu!**********
3325  Fmt: IMAGE "Lmd,Times,Count,Accu in Pincer.",3X,S.DDDDE,2(3X,DDD),3X,D.DDE    !******
3330      RETURN Lmd
3335    END IF
3340    Rendp=Rend
3345    Count=0
3350    Curaccu=Initaccu/Accufac^Times ! Determine accuracy for Tloop.
3355    Dl=Curaccu*Lmd
3360  ! ------------------------------------------------------
3365    FOR I=-1 TO 1
3370      L(I)=Lmd+I*Dl
3375      U(I)=FNU(M,N,C,L(I),Rend)
3380    NEXT I
3385  ! ------------------------------------------------------
3390  Iloop: Minii=-1            ! Find I which gives smallest |U|.
3395    Umin=U(-1)
3400    Pincerflg=1
3405    FOR I=0 TO 1
3410      IF U(I)<Umin THEN
3415        Umin=U(I)
3420        Minii=I
3425      END IF
3430    NEXT I
3435    Lmd=L(Minii)                  ! New Lamda.
3440    IF Minii=0 THEN
3445      Pincerflg=0
3450    ELSE
3455      Pincerflg=1
3460    END IF
3465    Count=Count+1                 ! Counter of Iloop.
3470    IF Curaccu<Accu OR (U(-1)=U(0) AND U(0)=U(1)) THEN
3475                      ! 2nd condition is needed because of EC accu.
3480      Times=Times+1               ! Counter of Tloop.
3485      GOTO Tloop
3490    END IF
3495  ! ------------------------------------------------------
3500  ! PRINT "FLG,MINII,DL",Pincerflg,Minii,Dl!******
3505    SELECT Pincerflg
3510    CASE 0                        ! When held between, increase accuracy.
3515      Curaccu=Curaccu/2
3520      Dl=Dl/2
3525      L(-1)=Lmd-Dl
3530      L(1)=Lmd+Dl
3535      L(0)=Lmd
3540      U(0)=U(Minii)
3545      FOR I=-1 TO 1 STEP 2
3550        U(I)=FNU(M,N,C,L(I),Rend)
3555      NEXT I
3560    CASE 1                        ! When not held between, shift.
3565      SELECT Minii
3570      CASE 1
3575        L(-1)=L(0)
3580        L(0)=L(1)
3585        L(1)=L(1)+Dl
3590        U(-1)=U(0)
3595        U(0)=U(1)
3600        U(1)=FNU(M,N,C,L(1),Rend)
3605      CASE -1
3610        L(1)=L(0)
3615        L(0)=L(-1)
3620        L(-1)=L(-1)-Dl
3625        U(1)=U(0)
3630        U(0)=U(-1)
3635        U(-1)=FNU(M,N,C,L(-1),Rend)
3640      END SELECT
3645    END SELECT
3650 ! PRINT USING "3(D.20DE,X)";U(1),U(0),U(-1)  !*********8
3655    DISP "Count of Iloop & Minii in Pincer ",Count,Minii
3660    GOTO Iloop
3665  FNEND
3670  ! ####################################################
3675  DEF FNU(M,N,C,L,Rend)
3680  ! #####################(D)############################
3685    RETURN ABS(FNNmrl(M,N-M+2,C,L)-FNNmru(M,N-M+2,C,L,Rend))
3690  FNEND
3695  ! ####################################################
3700  SUB Kappa(M,N,C,K,Kap1,Kap2)
3705  !  Joining factors. See Flammer eqs.(4.2.2)&(4.2.5).
3710  !    K     Kap1     Kap2
3715  !    1   Kappa(1)    0
3720  !    2      0     Kappa(2)
3725  !    3   Kappa(1) Kappa(2)
3730  !  d(-r) must be given for calculation of Kappa(2).
3735  ! #####################(C)############################
3740    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
3745  ! ------------------Summation---------------------------
3750    M2=2*M
3755    Mr=FNFact(M2+Rstart)
3760    S=Dd(Rstart)*Mr
3765    FOR R=Rstart+2 TO Rmax STEP 2
3770      Mr=Mr*(M2+R)*(M2+R-1)/R/(R-1)
3775      S=S+Dd(R)*Mr
3780    NEXT R
3785  ! ------------------------------------------------------
3790    Nmf=FNFact(N+M+Rstart)
3795    Nmmf=FNFact((N-M-Rstart)/2)
3800    Nmpf=FNFact((N+M+Rstart)/2)
3805    Mf=FNFact(M)
3810    SELECT K
3815    CASE 1
3820      Kap1=(M2+1+Rstart*2)*Nmf/Mf*S/2^(N+M)/Dd(Rstart)/C^(M+Rstart)/Nmmf/Nmpf
3825    CASE 2
3830      Kap2=(1-2*Rstart)*2^(N-M)*FNFact(M2)/Nmf*Nmmf*Nmpf/C^(M-1-Rstart)*Dd(-M2+Rstart)*S/(1+Rstart*(M2-4))/(M2-1)/Mf
3835    CASE 3
3840      Kap1=(M2+1+Rstart*2)*Nmf/Mf*S/2^(N+M)/Dd(Rstart)/C^(M+Rstart)/Nmmf/Nmpf
3845      Kap2=(1-2*Rstart)*2^(N-M)*FNFact(M2)/Nmf*Nmmf*Nmpf/C^(M-1-Rstart)*Dd(-M2+Rstart)*S/(1+Rstart*(M2-4))/(M2-1)/Mf
3850    END SELECT
3855  SUBEND
3860  ! ####################################################
3865  DEF FNLamdalc(M,N,C)
3870  !  Lamda for large C.  See Flammer eq.(8.1.11).
3875  ! #####################(L)############################
3880    READ C1,C2                 ! Interporation bound for m=0,n=3.
3885    DATA 10,10.15
3890    IF M=0 AND N=3 AND C>C1 AND C<C2 THEN  ! Interporation.
3895      Lc1=62.2577002498796
3900      Lc2=63.3325655538975
3905      RETURN Lc1+(Lc2-Lc1)*(C-C1)/(C2-C1)
3910    END IF
3915    DIM L(6),L2(7),Mm(6),T(-20:16),Cc(-5:0)
3920    Ll=N-M
3925    L(1)=Ll
3930    FOR I=2 TO 6
3935      L(I)=L(I-1)*Ll            ! L(I)=L^I
3940    NEXT I
3945    L2(1)=2*Ll+1
3950    FOR I=3 TO 7 STEP 2
3955      L2(I)=L2(I-2)*L2(1)*L2(1) ! L2(I)=(2L+1)^I
3960    NEXT I
3965    Mm(2)=M*M
3970    FOR I=4 TO 6 STEP 2
3975      Mm(I)=Mm(I-2)*Mm(2)       ! Mm(I)=M^I
3980    NEXT I
3985    T(1)=2
3990    FOR I=2 TO 16
3995      T(I)=T(I-1)*2             ! T(I)=2^I
4000    NEXT I
4005    T(0)=1
4010    FOR I=-1 TO -20 STEP -1
4015      T(I)=T(I+1)/2
4020    NEXT I
4025    Cc(0)=1
4030    FOR I=-1 TO -5 STEP -1
4035      Cc(I)=Cc(I+1)/C           ! Cc(I)=C^I
4040    NEXT I
4045  ! ------------------------------------------------------
4050    A=L2(1)*C-(2*L(2)+2*Ll+3-4*Mm(2))*T(-2)-L2(1)*(L(2)+Ll+3-8*Mm(2))*T(-4)*Cc(-1)
4055    B=-T(-6)*Cc(-2)*(5*(L(4)+2*L(3)+7*Ll+3)-48*Mm(2)*(2*L(2)+2*Ll+1))
4060    D=66*L(5)+165*L(4)+962*L(3)+1278*L(2)+1321*Ll+453
4065    E=-Mm(2)*(2368*L(3)+3552*L(2)+4448*Ll+1632)+Mm(4)*(256*Ll+128)
4070    F=252*L(6)+756*L(5)+5885*L(4)+10510*L(3)+18478*L(2)+13349*Ll+4425
4075    G=-Mm(2)*(14720*L(4)+29440*L(3)+64000*L(2)+49280*Ll+17280)+Mm(4)*(6144*L(2)+6144*Ll+3072)
4080    H=527*L2(7)+61529*L2(5)+1043961*L2(3)+2241599*L2(1)
4085    J=T(5)*Mm(2)*(5739*L2(5)+127550*L2(3)+298951*L2(1))
4090    K=-T(11)*Mm(4)*(355*L2(3)+1505*L2(1))+T(16)*Mm(6)*L2(1)
4095    Lamda=A+B-(D+E)*T(-10)*Cc(-3)-(F+G)*T(-12)*Cc(-4)-(H+J+K)*T(-20)*Cc(-5)
4100    ! ------------------Correction--------------------------
4105    IF M=0 THEN
4110      IF N=3 AND C>10.18 AND C<=10.3 THEN Lamda=Lamda+.1
4115        ! 10<C<10.15 Lamda is not correct.
4120      IF N=4 AND C>9.6 AND C<10 THEN Lamda=Lamda-.25
4125      IF N=5 AND C>10.1 AND C<10.4 THEN Lamda=Lamda-.6
4130      IF N=6 AND C>10.7 AND C<11.1 THEN Lamda=Lamda-.73
4135      IF N=7 AND C>11.5 THEN Lamda=Lamda-.5
4140    END IF
4145    RETURN Lamda
4150  FNEND
4155  ! ####################################################
4160  DEF FNLamdasc(M,N,C)
4165  !  Lamda by power series of c. See Flammer eq.(3.1.17).
4170  !  Use FNL.  Lpre must be called beforehand.
4175  ! #####################(L)############################
4180    COM /Lv/Lv(*)            ! Common with Lpre.
4185    L=0
4190    IF C=0 THEN C=1.E-10
4195    FOR K2=0 TO 10 STEP 2
4200      L=L+Lv(K2)*C^K2
4205    NEXT K2
4210    RETURN L
4215  FNEND
4220  ! ####################################################
4225  DEF FNNmru(M,R,C,L,Rend)
4230  !   Nm,r when r ascends.  See Flammer eqs.(3.1.8) & (3.1.11).
4235  !   U2 is -Nmru whose R is replaced by n-m+2. See Flammer eq.(3.1.15).
4240  ! #####################(D)############################
4245    COM /C/C4
4250    SELECT R
4255    CASE >Rend
4260      RETURN 0                    ! Stop recursion.
4265    CASE <=Rend
4270    ! Pgamma(M,R,C,Pg)
4275    ! Pbc(M,R,Pbc)
4280    ! C2=C*C
4285    ! RETURN -Pbc*C4/(Pg-L+FNNmru(M,R+2,C,L,Rend))
4290      RETURN -FNBc(M,R)*C4/(FNGamma(M,R,C)-L+FNNmru(M,R+2,C,L,Rend))
4295    END SELECT
4300  FNEND
4305  ! ####################################################
4310  SUB Rendfind(M,R,C,L,Rend)
4315  ! #####################(D)############################
4320    DATA 1E-10,4
4325    READ Errormin,Rstep             ! Errormin changes accu.
4330    Error=1
4335    Rend=R+2
4340    WHILE Error>Errormin
4345      A=FNNmru(M,R,C,L,Rend)
4350      B=FNNmru(M,R,C,L,Rend+2)
4355      Error=ABS((A-B)/B)
4360      Rend=Rend+Rstep
4365    END WHILE
4370  SUBEND
4375  ! ####################################################
4380  DEF FNNmrl(M,Rp2,C,L)
4385  !   Nm,r when r descends.  See Flammer eqs.(3.1.9) & (3.1.12).
4390  !   U1 is Nmrl whose Rp2 is replaced by n-m+2.  See Flammer eq.(3.1.14).
4395  ! #####################(D)############################
4400    COM /C/C4
4405    R=Rp2-2
4410    SELECT R
4415    CASE -.1 TO 1.1
4420  !   Pgamma(M,R,C,Pg)
4425  !   RETURN -Pg+L
4430      RETURN -FNGamma(M,R,C)+L        ! Stop recursion.
4435    CASE >1.9
4440  !   Pgamma(M,R,C,Pg)
4445  !   Pbc(M,R,Pbc)
4450  !   RETURN -Pg+L-Pbc*C4/FNNmrl(M,R,C,L)
4455      RETURN -FNGamma(M,R,C)+L-FNBc(M,R)*C4/FNNmrl(M,R,C,L)
4460    END SELECT
4465  FNEND
4470  ! ####################################################
4475  DEF FNLgna1(M,N,R,X)
4480  !  Asociated Legendre function of first kind.
4485  !  See "Handbook of value analysis by FORTRAN" by Isoda & Ohno,
4490  !   Ohm Co.,Tokyo,1971,p508.
4495  !    R=0 when X real, 1 imaginary.
4500  !    If M=0 then Legendre func. of 1st kind.
4505  ! #######################(SF-2)#######################
4510    INTEGER I,J,K,Nm,Nn,Mn,Mi,N1,N2,N3
4515    DIM G(1:200)                  ! Factorial. Must be 10^308
4520    IF N=0 AND M=0 THEN RETURN 1  ! P0,0(x)=1.
4525    Nm=N-M
4530    Z=1
4535    W=Z
4540    SELECT Nm
4545    CASE <0
4550      RETURN 0                    ! Pm,0(x)=0 (m<>0)
4555    CASE >0
4560      FOR I=1 TO Nm
4565        Z=X*Z
4566        IF ABS(Z)<1.E-200 THEN
4567          Z=0
4568          GOTO 4575
4569        END IF
4571      NEXT I
4575    END SELECT
4580    G(1)=1
4585    Nn=N+N+1
4590    FOR I=2 TO Nn
4595      G(I)=W*G(I-1)
4600      W=W+1
4605    NEXT I
4610    W=1
4615    Y=W/(X*X)
4620    IF R=0 THEN
4625      Y=-Y
4630      W=-W
4635    END IF
4640    SELECT X
4645    CASE 0
4650      I=Nm/2
4655      IF 2*I-Nm<>0 THEN RETURN 0
4660      Mn=M+N+1
4665      Mi=M+I+1
4670      P=G(Mn)/(G(I+1)*G(Mi))
4675      K=I/2
4680      IF 2*K-I<>0 THEN P=-P
4685    CASE <>0
4690      J=3
4695      P=0
4700      FOR I=1 TO 12
4705        IF (Nm+2)/2-I<0 THEN GOTO L400
4710        N1=N+N-I-I+3
4715        N2=N-I+2
4720        N3=N-I-I-M+J
4725        P=P+G(N1)*Z/(G(I)*G(N2)*G(N3))
4730        Z=Z*Y
4735      NEXT I
4740    END SELECT
4745  L400: Z=1
4750    FOR I=1 TO N
4755      Z=Z+Z
4760    NEXT I
4765    P=P/Z
4770    IF R=0 THEN GOTO L520
4775    I=N-N/4
4780    IF I-1>0 THEN P=-P
4785  L520: IF M=0 THEN RETURN P
4790    J=M/2
4795    Z=ABS(W+X*X)
4800    IF M-2*J=0 THEN GOTO L610
4805    Z=SQR(Z)
4810    J=M
4815  L610: FOR I=1 TO J
4820      P=P*Z
4825    NEXT I
4830    RETURN P
4835  FNEND
4840  ! ####################################################
4845  SUB Lgna2(Mmax,Nmax,Nmin,X,Q(*))
4850  !  Associated Legendre function of the 2nd kind.
4855  !  See "Handbook of value analysis" by Isoda & Ohno, Ohm Co.,1971,Tokyo.
4860  !    Caution! There are some errors.
4865  !  Results are stored in matrix with sufixes of M,N.
4870  !  For only +n, Nmin=0.
4875  !  For -n, Nmin=-n & Nmax>=1.
4880  !  Qmn is infinite when n<-m.
4885  ! ####################################################
4890    ! ------------------Error-------------------------------
4895    IF X<=1 THEN
4900      BEEP
4905      PRINT "X<=1 in Lgna2. Return 0."
4910      SUBEXIT
4915    END IF
4920    IF Nmax<0 THEN
4925      BEEP
4930      PRINT "Nmax<0 in Lgna2. Return 0."
4935      SUBEXIT
4940    END IF
4945    IF Mmax<0 THEN
4950      BEEP
4955      PRINT " Mmax<0 in Lgna2. Return 0."
4960      SUBEXIT
4965    END IF
4970    IF Nmin>0 THEN
4975      BEEP
4980      PRINT "Nmin>0 in Lgna2.  Nmin changed to 0."
4985      Nmin=0
4990    END IF
4995    IF Nmin<0 AND Nmax<1 THEN
5000      BEEP
5005      PRINT "Nmin<0 & Nmax<1 in Lgna2.  Nmax changed to 1."
5010      Nmax=1
5015    END IF
5020    REDIM Q(Mmax,Nmin:Nmax)
5025    ! ------------------Initial value-----------------------
5030    Tau=.5*LOG(1+2/(X-1))
5035    Q(0,0)=Tau
5040    IF Nmax>0 THEN Q(0,1)=X*Tau-1
5045    Rx=SQR(X*X-1)
5050    IF Mmax>0 THEN Q(1,0)=-1/Rx
5055    IF Mmax*Nmax>=1 THEN Q(1,1)=(X*Q(0,1)-Q(0,0))/Rx  ! Eq.(8.8.11)
5060    ! ----------------Recursion in n------------------------
5065    IF Nmax>=2 THEN
5070      FOR M=0 TO 1*(Mmax<>0)
5075        FOR N=2 TO Nmax
5080          Q(M,N)=((2*N-1)*X*Q(M,N-1)-(N+M-1)*Q(M,N-2))/(N-M)
5085        NEXT N
5090      NEXT M
5095    END IF
5100    ! ----------------Recursion in m------------------------
5105    IF Mmax>=2 THEN
5110      FOR N=0 TO Nmax
5115        FOR M=2 TO Mmax
5120          Q(M,N)=(N-M+2)*(N+M-1)*Q(M-2,N)-2*(M-1)*X/Rx*Q(M-1,N)
5125        NEXT M
5130      NEXT N
5135    END IF
5140    IF Nmin=0 THEN SUBEXIT      !  No negative n.
5145    ! ------------------Negative n--------------------------
5150    !  Backward recursion in n.
5155    Inf=1.E+20
5160    FOR M=0 TO Mmax
5165      FOR N=-1 TO Nmin STEP -1
5170        IF N<-M THEN
5175          Q(M,N)=Inf
5180        ELSE
5185          Q(M,N)=((2*N+3)*X*Q(M,N+1)-(N-M+2)*Q(M,N+2))/(N+M+1)
5190        END IF
5195      NEXT N
5200    NEXT M
5205  SUBEND
5210  ! ####################################################
5215  DEF FNSnbes(Nn,Xx)
5220    !  Spherical Bessel function of the second kind.
5225    !  See "Value Analysis and FORTRAN" by Amamiya & Taguchi,
5230    !    Maruzen Co.,Tokyo,1971.
5235    ! #####################(SF-6)#########################
5240    N=Nn
5245    X=Xx
5250    IF N>=3.E+4 THEN GOSUB Notaccu
5255    IF X<>0 THEN Z=1/X
5260    SELECT X
5265    CASE <0
5270      GOSUB Invalid
5275    CASE 0
5280      RETURN -1.E+30
5285    CASE 0 TO 7.E-4
5290      GOSUB Small
5295    CASE ELSE
5300      GOSUB Large
5305    END SELECT
5310    ! ------------------------------------------------------
5315  Notaccu: PRINT "Value of SNBES is not accurate.  N=";N;"  X=";X
5320    RETURN 0
5325  Invalid: PRINT "Argument of SNBES is invalid.  N=";N;"  X=";X
5330    RETURN 0
5335  Over: PRINT "Value of SNBES is overflow.  N=";N;"  X=";X
5340    RETURN -1.E+30
5345    ! ------------------------------------------------------
5350  Small: SELECT N
5355    CASE <0
5360      GOSUB Invalid
5365    CASE 0
5370      RETURN -Z
5375    CASE >0
5380      M=30/LGT(Z)
5385      IF N-M+1>0 THEN GOSUB Over
5390      Qn0=Z
5395      Qn1=1
5400      FOR I=1 TO N
5405        Qn2=Qn0*Qn1*Z
5410        Qn1=Qn1+2
5415        Qn0=Qn2
5420      NEXT I
5425      RETURN -Qn2
5430    END SELECT
5435    ! ------------------------------------------------------
5440  Large: Qn0=-Z*COS(X)
5445    Qn1=Z*(Qn0-SIN(X))
5450    SELECT N
5455    CASE <0
5460      GOSUB Invalid
5465    CASE 0
5470      RETURN Qn0
5475    CASE 1
5480      RETURN Qn1
5485    CASE >1
5490      FOR I=2 TO N
5495        Qn2=(I+I-1)*Z*Qn1-Qn0
5500        IF Qn2+1.E+30<=0 THEN GOSUB Over
5505        Qn0=Qn1
5510        Qn1=Qn2
5515      NEXT I
5520      RETURN Qn2
5525    END SELECT
5530  FNEND
5535  ! ####################################################
5540  DEF FNFact(N)
5545  !  Factorial of N, N!.
5550  ! #######################(SF-8)#######################
5555    SELECT N
5560    CASE <0
5565      PRINT "Invalid data in FNFactorial. N<0. Return 1."
5570      RETURN 1
5575    CASE 0
5580      RETURN 1
5585    CASE ELSE
5590      F=1
5595      FOR I=N TO 1 STEP -1
5600        F=F*I
5605      NEXT I
5610      RETURN F
5615    END SELECT
5620  FNEND
5625  ! ####################################################
5630  SUB Sbes_der(N,X,Sb$,K,Sb0,Sb1,Sb2)
5635  ! 1st & 2nd derivatives of spherical Bessel functions.
5640  ! See Abramowitz & Stegun eq.10.1.24.
5645  !   Sb$      J:s_Bessel, N:s_Neumann
5650  !   K        1        2
5655  !   Sb0      Sb       Sb
5660  !   Sb1      Sb'      Sb'
5665  !   Sb2      0        Sb''
5670  ! Need FNSjbes or Snbes.
5675  ! #####################(SF-9)#########################
5680  !DISP "Now in Sbes_der."
5685  ! ------------------1st derivative----------------------
5690    Sb2=0
5695    SELECT Sb$
5700    CASE "J"
5705      Sb0=FNSjbes(N,X)
5710      Sb01=FNSjbes(N+1,X)
5715    CASE "N"
5720      Sb0=FNSnbes(N,X)
5725      Sb01=FNSnbes(N+1,X)
5730    END SELECT
5735    Sb1=N/X*Sb0-Sb01
5740    IF K=1 THEN SUBEXIT
5745  ! ------------------2nd derivative----------------------
5750    IF Sb$="J" THEN Sb02=FNSjbes(N+2,X)
5755    IF Sb$="N" THEN Sb02=FNSnbes(N+2,X)
5760    Sb2=(2*N+1)/X*Sb1-N*(N+2)/X/X*Sb0+Sb02
5765  SUBEND
5770  ! ####################################################
5775  SUB Lgna_der(M,N,X,Pq$,Lgna,Lgnad)
5780  !  Derivative of associated Legendre functions.
5785  !  Pq$      P: 1st kind,   Q: 2nd kind
5790  !  Lgna     Associated Legendre function (result).
5795  !  Lgnad    Its derivative (result).
5800  !  See Isoda & Ohno eq.(8.8.2).
5805  ! #####################(SF-10)########################
5810    SELECT Pq$
5815    CASE "P"
5820      Lgna=FNLgna1(M,N,0,X)
5825      Lgna1=FNLgna1(M,N+1,0,X)
5830    CASE "Q"
5835      SELECT N
5840      CASE <0
5845        Nmax=1                      ! Use 0 & 1 order for negative N.
5850        Nmin=N
5855      CASE >=0
5860        Nmax=N+1                    ! Use N+1 order for derivative.
5865        Nmin=0
5870      END SELECT
5875      Lgna2(M,Nmax,Nmin,X,Q(*))
5880      Lgna=Q(M,N)
5885      Lgna1=Q(M,N+1)
5890    END SELECT
5895  ! PRINT "Lgna1,Lgna in Lgna_der",Lgna1,Lgna
5900    Lgnad=((M-N-1)*Lgna1+(N+1)*X*Lgna)/(1-X*X)
5905  SUBEND
5910  ! ####################################################
5915  DEF FNCabs(R,I)
5920  ! ####################################################
5925    RETURN SQR(R*R+I*I)
5930  FNEND
5935  ! ####################################################
5940  SUB Cdivid(Xr,Xi,Yr,Yi,Rr,Ri)
5945  ! Division of complex number.
5950  ! Rr+iRi=(Xr+iXi)/(Yr+iYi)
5955  ! ####################################################
5960    Inf=1.E+30
5965    C=Yr*Yr+Yi*Yi
5970    SELECT C
5975    CASE 0
5980      BEEP
5985      DISP "Error in Cdivid. Divisor is 0. Return ";Inf
5990      WAIT 1
5995      Rr=0
6000      Ri=0
6005      SUBEXIT
6010    CASE ELSE
6015      Rr=(Yr*Xr+Yi*Xi)/C
6020      Ri=(Yr*Xi-Yi*Xr)/C
6025    END SELECT
6030  SUBEND
6035  ! ####################################################
6040  SUB Cinv(A(*),N1,N2,Det(*),Ind)
6045  ! Inversion of complex matrix and solution of simultaneous equations.
6050  ! See Amamiya & Taguchi "Value Analysis and FORTRAN",Maruzen,s44.
6055  ! Revised for complex number by Furusawa.
6060  ! A(I,J,K)   N1*(N1+N2)*2 dimension array.
6065  ! N1         Size of the (coefficient) array.
6070  ! N2         N2 constants vectors for simultaneouse equations.
6075  !            N2=0 for only inversion, N2=1 for one set of constants.
6080  ! Det(K)     Determinant of A.
6085  ! Ind        Status indicator.
6090  !            1=(pivot>10^-4), 2=(10^-6 to 10^-4), 3=(<10^-6) and error
6095  !            4=imput error
6100  ! ####################################################
6105    OPTION BASE 1
6110    DIM Pivot(2),Pivi(2),Perm(50),X(100,2)
6115    READ Singval
6120    DATA 1E+38
6125    N=N1
6130    M=N+N2
6135    IF N<=0 OR N2<0 OR N>50 OR M>100 THEN
6140      BEEP
6145      PRINT "Input error in Cinv.  0<N1<=50, 0<=N2 and N1+N2<=100 must be satisfied. Returns with no calculation and Ind=4."
6150      Ind=4
6155      SUBEXIT
6160    END IF
6165    Ind=1
6170    Det(1)=1
6175    Det(2)=0
6180    FOR I=1 TO N
6185      Perm(I)=I                ! Memory for re-arrangement of column.
6190    NEXT I
6195    Eps=0
6200    ! ------------------------------------------------------
6205    FOR K=1 TO N
6210      Max=0                    ! Pivot is selected for max. value of column.
6215      FOR J=K TO N
6220        W=FNCabs(A(K,J,1),A(K,J,2))
6225        IF Max<=W THEN
6230          Max=W
6235          L=J                  ! L is pivot's column.
6240        END IF
6245      NEXT J
6250      SELECT Max
6255      CASE >Eps                ! OK
6260        GOTO Pivot
6265      CASE Eps*.01 TO Eps      ! Warning
6270        Ind=2
6275        GOTO Pivot
6280      CASE ELSE                ! Singular
6285        MAT Det=(0)
6290        Ind=3
6295        MAT A=(1.E+38)
6300        BEEP
6305        PRINT "Given matrix is singular in Cinv. Returned values =";Singval
6310        SUBEXIT
6315      END SELECT
6320      ! ------------------------------------------------------
6325  Pivot: Pivot(1)=A(K,L,1)
6330      Pivot(2)=A(K,L,2)
6335      Cmult(Det(1),Det(2),Pivot(1),Pivot(2),Detr(1),Detr(2))
6340      Det(1)=Detr(1)
6345      Det(2)=Detr(2)
6350      Cdivid(1,0,Pivot(1),Pivot(2),Pivi(1),Pivi(2))
6355      IF L=K THEN Sweep
6360      Iw=Perm(K)               ! Replacement between K and L.
6365      Perm(K)=Perm(L)
6370      Perm(L)=Iw
6375      FOR I=1 TO N
6380        Creplace(A(I,K,1),A(I,K,2),A(I,L,1),A(I,L,2))
6385      NEXT I
6390      ! ------------------------------------------------------
6395  Sweep: FOR J=1 TO M
6400        Cmult(A(K,J,1),A(K,J,2),Pivi(1),Pivi(2),X(J,1),X(J,2))
6405        A(K,J,1)=X(J,1)
6410        A(K,J,2)=X(J,2)
6415      NEXT J
6420      FOR I=1 TO N
6425        IF I=K THEN Nexti
6430        W1=A(I,K,1)
6435        W2=A(I,K,2)
6440        IF W1=0 AND W2=0 THEN Nexti
6445        FOR J=1 TO M
6450          IF J=K THEN Nextj
6455          Cmult(W1,W2,X(J,1),X(J,2),W11,W22)
6460          A(I,J,1)=-W11+A(I,J,1)
6465          A(I,J,2)=-W22+A(I,J,2)
6470  Nextj: NEXT J
6475        Cmult(-W1,-W2,Pivi(1),Pivi(2),A(I,K,1),A(I,K,2))
6480  Nexti: NEXT I
6485      A(K,K,1)=Pivi(1)
6490      A(K,K,2)=Pivi(2)
6495      Eps=MAX(Max*1.E-4,Eps)
6500    NEXT K
6505    ! ------------------------------------------------------
6510    FOR I=1 TO N
6515  Again: K=Perm(I)
6520      IF K=I THEN Skip
6525      Iw=Perm(K)
6530      Perm(K)=Perm(I)
6535      Perm(I)=Iw
6540      FOR J=1 TO M
6545        Creplace(A(I,J,1),A(I,J,2),A(K,J,1),A(K,J,2))
6550      NEXT J
6555      Det(1)=-Det(1)
6560      Det(2)=-Det(2)
6565      GOTO Again
6570  Skip: NEXT I
6575  SUBEND
6580  ! ####################################################
6585  SUB Cmult(Xr,Xi,Yr,Yi,Rr,Ri)
6590  ! ####################################################
6595    Rr=Xr*Yr-Xi*Yi
6600    Ri=Xr*Yi+Xi*Yr
6605  SUBEND
6610  ! ####################################################
6615  SUB Creplace(Xr,Xi,Yr,Yi)
6620  ! Replacement of complex numbers Xr + i Xi and Yr + i Yi.
6625  ! ####################################################
6630    Wr=Xr
6635    Wi=Xi
6640    Xr=Yr
6645    Xi=Yi
6650    Yr=Wr
6655    Yi=Wi
6660  SUBEND
6665  ! ####################################################
6670  SUB Label_length(X$,Xn)
6675  ! Check label length. LEN(X$) must be <=Xn.
6680  ! #######################(GP-8)#######################
6685    Xnc=LEN(X$)
6690    IF Xnc>Xn THEN
6695      BEEP
6700      PRINT "Current charactor number is ";Xnc;".  Shorten by ";Xnc-Xn;" charactors."
6705      DISP "Current label: ",X$
6710      INPUT "Shotened label ?",X$
6715    END IF
6720  SUBEND
6725   ! ####################################################
6730  Label: SUB Label(Text$)
6735   ! ######################(GP-9)########################
6740   !  This prints a character string at the current pen position and using
6745   !  the current LORG, LDIR and CSIZE.  The LORG will need to be redeclared
6750   !  upon returning to the calling context, as this routine needs LORG 1 if
6755   !  the text is longer than one character.
6760    OPTION BASE 1
6765    COM /Udc/Old_chars$[50],Size(50),Chars(50,40,3)
6770    REAL Array(41,3)
6775    FOR Char=1 TO LEN(Text$)
6780      IF Char=2 THEN LORG 1  ! Necessary when doing one character at a time
6785      Char$=Text$[Char;1]
6790      Pos=POS(Old_chars$,Char$)           ! Is this to be replaced by a UDC?
6795      IF Pos THEN
6800        REDIM Array(Size(Pos),3)                   ! \
6805        FOR Row=1 TO Size(Pos)                     !  \   Take a slice out
6810          FOR Column=1 TO 3                        !   >  of the 3D array
6815            Array(Row,Column)=Chars(Pos,Row,Column)!  /   and put it in the
6820          NEXT Column                              ! /    2D array for
6825        NEXT Row                                   !/     SYMBOL.
6830        WHERE X,Y
6835        SYMBOL Array(*)
6840        MOVE X,Y
6845        LABEL USING "#,K";" " ! Tell the computer to update the pen position
6850      ELSE ! (regular character)
6855        LABEL USING "#,K";Char$
6860      END IF ! (this character been redefined?)
6865    NEXT Char
6870  SUBEND
6875   ! ####################################################
6880  SUB Arcdeg(X,Y,P,D,R,Sd,Wd)          ! Radius,StartDeg,WidthDeg
6885  ! ####################################################
6890    MOVE X,Y                         ! Move in current PIVOT
6895    PIVOT D+Sd                       ! Total PIVOT
6900    PEN P
6905    POLYLINE R,360,Wd
6910    PIVOT D                          ! Return to original PIVOT
6915  SUBEND
6920  ! ####################################################
6925  SUB Densf
6930  ! Normalized TS of liquid versus Rho1by0 with parameter C1/0.
6935  ! Can compare Scatfunl(exac.) & Scatfunslr(appro.).
6940  ! ####################################################
6945    INPUT "Method(Appro./Exac.) ?",Ae$
6950    PRINT "Method(Appro./Exac.)",Ae$
6955    INPUT "Mmax,Nmax ?",Mmax,Nmax
6960    PRINT "Mmax,Nmax",Mmax,Nmax
6965    IF Ae$="A" THEN
6970      DIM Sf(2,0,1)                ! Kind,Theta,Phi
6975    ELSE
6980      DIM Sfe(0,1)                  ! Theta,Phi
6985    END IF
6990    INPUT "K0a, b/a ?",K0a,Bbya
6995    PRINT "K0a, b/a",K0a,Bbya
7000    INPUT "Rho1by0: Min,Max,Step ?",Rhomin,Rhomax,Rhostep
7005    PRINT "Rho1by0: Min,Max,Step",Rhomin,Rhomax,Rhostep
7010    INPUT "C1/0(para.): Min,Max,Step ?",Cmin,Cmax,Cstep
7015    PRINT "C1/0(para.): Min,Max,Step",Cmin,Cmax,Cstep
7020    ! ------------------Pre cal.----------------------------
7025    Xi0=1/SQR(1-Bbya^2)
7030    H0=K0a*Xi0
7035    IF Cstep=0 THEN Cstep=1
7040    IF Rhostep=0 THEN Rhostep=1
7045    Cimax=INT((Cmax-Cmin)/Cstep+.1)
7050    Rimax=INT((Rhomax-Rhomin)/Rhostep+.1)
7055    ALLOCATE F(Cimax,Rimax),Linlbl$(Cimax)[10]
7060    Prelfile(Paraok)
7065    T1=TIMEDATE
7070  ! ------------------------------------------------------
7075    Ci=0
7080    FOR C=Cmin TO Cmax STEP Cstep
7085      K1a=K0a/C
7090      H1=K1a*Xi0
7095      Ri=0
7100      FOR Rho1by0=Rhomin TO Rhomax STEP Rhostep
7105        IF Ae$="A" THEN
7110          Scatfunslr(1,Mmax,Nmax,H0,H1,Rho1by0,Xi0,90,-1,0,180,Sf(*))
7115          IF Sf(1,0,0)<=0 THEN Sf(1,0,0)=1.E-20
7120          F(Ci,Ri)=20*LGT(Sf(1,0,0)/K0a/2)
7125        ELSE
7130          Scatfunl(Mmax,Nmax,H0,H1,Rho1by0,Xi0,90,-1,0,180,Sfe(*))
7135          IF Sfe(0,0)<=0 THEN Sfe(0,0)=1.E-20
7140          F(Ci,Ri)=20*LGT(Sfe(0,0)/K0a/2)
7145        END IF
7150        BEEP
7155        PRINT "C1/0 ,Rho1by0, Norm.TS, Time",C,Rho1by0,DROUND(F(Ci,Ri),3),DROUND((TIMEDATE-T1)/3600,2)
7160        Ri=Ri+1
7165      NEXT Rho1by0
7170      Linlbl$(Ci)=VAL$(C)
7175      Ci=Ci+1
7180    NEXT C
7185    T2=TIMEDATE
7190    BEEP 2000,3
7195    PRINT "Time",T2-T1
7200  ! ------------------Plot--------------------------------
7205    DIM Led$[80]
7210    Led$="DF "&Ae$&" M="&VAL$(Mmax)&" N="&VAL$(Nmax)&" b/a="&VAL$(Bbya)
7215    Led$=Led$&" Xi0="&VAL$(DROUND(Xi0,5))&" k0a="&VAL$(DROUND(K0a,3))
7220    Plot(F(*),Rhomin,Rhomax,Rhostep,"LIN","Rho1/0","20log(g/2a)"," ",0,Led$,Linlbl$(*),1)
7225    GOTO 6945
7230  SUBEND
7235  ! ####################################################
7240  SUB Rmnhc
7245     ! Rmnh check. Comparison with Hanish et al.
7250     ! ####################################################
7255    COM /Pltr/Plot$[8],Sp$[8],Lastpen
7260    DIM Led$[80]
7265    Prelfile(Paraok)
7270    INPUT "Graph(G) or Value(V) ?",Gv$
7275    SELECT Gv$
7280    CASE "G"
7285  In1: INPUT "M,N,C",M,N,Cc
7290  Ka: INPUT "Kamin,Kamax,Kainc",Kamin,Kamax,Kainc
7295      IF Kamin<Cc THEN
7300        PRINT "Kamin must be lager than C. C=";Cc
7305        GOTO Ka
7310      END IF
7315      Imax=INT((Kamax-Kamin)/Kainc)
7320      ALLOCATE R(1:4,0:Imax)
7325      I=0
7330      FOR Ka=Kamin TO Kamax STEP Kainc
7335        DISP "Ka",Ka
7340        Presphe(M,N,Cc,L)
7345        Xi=Ka/Cc
7350        Rmnh(Kind,M,N,Cc,L,Xi,R(1,I),R(2,I),R(3,I),R(4,I),Err)
7355        I=I+1
7360      NEXT Ka
7365      BEEP 2000,3
7370        ! ------------------------------------------------------
7375      DIM Linlbl$(0)[10]
7380      Led$="Rmnhc  "&" M="&VAL$(M)&" N="&VAL$(N)&" C="&VAL$(Cc)
7385  Plot: Plot(R(*),Kamin,Kamax,Kainc,"LIN","ka","R","Hanish",0,Led$,Linlbl$(*),1)
7390      INPUT "Plot or Other case (P/C)",Pc$
7395      IF Pc$="P" THEN GOTO Plot
7400      IF Pc$="C" THEN
7405        DEALLOCATE R(*)
7410        GOTO In1
7415      END IF
7420    CASE "V"
7425      DIM C(10)
7430      INPUT "Printer ?(P/C)",Pr$
7435      IF Pr$="P" THEN
7440        PRINTER IS PRT
7445      ELSE
7450        PRINTER IS CRT
7455      END IF
7460  In: INPUT "Xi,Mmin,Mmax, Nmax ?",Xi,Mmin,Mmax,Nmax
7465      I=0
7470      INPUT "C(0 for stop)",C(I)
7475      I=I+1
7480      IF C(I-1)<>0 THEN GOTO 7470
7485      Imax=I-2
7490      FOR I=0 TO Imax
7495        PRINT
7500        PRINT "C=";C(I),"Xi=";Xi
7505        PRINT
7510        PRINT "M  N         R1        R1D         R2        R2D         L        ERROR"
7515        PRINT
7520        T1=TIMEDATE
7525        FOR M=Mmin TO Mmax
7530          FOR N=M TO Nmax
7535            Presphe(M,N,C(I),L)
7540     !       Dminus(M,N,C(I),L)
7545            Rmnh(30,M,N,C(I),L,Xi,R1,R2,R1d,R2d,Err)
7550            PRINT USING "2(DD,X),6(SD.3DE,1X)";M,N,R1,R1d,R2,R2d,L,Err
7555          NEXT N
7560        NEXT M
7565        T2=TIMEDATE
7570        PRINT "Time=";T2-T1
7575      NEXT I
7580      BEEP 1000,2
7585      GOTO In
7590    END SELECT
7595    PRINTER IS CRT
7600  SUBEND
7605  ! ####################################################
7610  SUB Rmn_c2
7615    ! Rmn check for fixed c & xi and varied m & n.
7620    ! ####################################################
7625    PRINTER IS PRT
7630    INPUT "Xi,Mmax, Nmax ?",Xi,Mmax,Nmax
7635    DIM C(10)
7640    I=0
7645    INPUT "C(0 for stop)",C(I)
7650    I=I+1
7655    IF C(I-1)<>0 THEN GOTO 7645
7660    Imax=I-2
7665    FOR I=0 TO Imax
7670      PRINT
7675      PRINT "C=";C(I),"Xi=";Xi
7680      PRINT
7685      PRINT "M  N      1JF      1SB      2JF      2SB      1JF'     1SB'     2JF'     2SB'"
7690      PRINT
7695      T1=TIMEDATE
7700      FOR M=0 TO Mmax
7705        FOR N=M TO Nmax
7710          Presphe(M,N,C(I),L)
7715          Dminus(M,N,C(I),L)
7720          Rmnjf(M,N,C(I),Xi,30,D,R1jf,R2jf,R1djf,R2djf)
7725          Rmnsb(M,N,C(I),Xi,30,R1sb,R2sb,R1dsb,R2dsb)
7730          PRINT USING "2(DD,X),8(SD.2DE)";M,N,R1jf,R1sb,R2jf,R2sb,R1djf,R1dsb,R2djf,R2dsb
7735        NEXT N
7740      NEXT M
7745      T2=TIMEDATE
7750      PRINT "Time=";T2-T1
7755    NEXT I
7760    BEEP 1000,2
7765    PRINTER IS CRT
7770  SUBEND
7775  ! ####################################################
7780  SUB Dbyd2(M,N,C,L)
7785  !  Revised from Dbyd.
7790  !  d(m,n,r)/d(m,n,n-m). See Flammer eq.(3.1.7)-(3.1.10) & 3.1.3.
7795  !  Needs FNLamda,(P)Gamma,(P)Bc and Nd.
7800  !  First, lamda is calculated with Laccu.
7805  !  Rmax is determined to give minimum d/d or truncated to Rmaxmax.
7810  !  FNLamda or FNLamdaf is selected by COM /Lcf/.
7815  !  L is used in this SUB and transfered to Dminus.
7820  ! #####################(D)############################
7825    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
7830                              ! Dd,Rstart & Rmax  are results.
7835    COM /Lcf/Lcf$
7840    COM /C/C4
7845    READ Laccu,Rmaxmax,Rmaxbias,Daccu  ! Daccu is only for caution.
7850    DATA 1E-15,50,5,0.01     ! Laccu must be determined by
7855                              !   Current_accu in Bouwkamp (EC accu).
7860    ! --------------Intial Rmax & L--------------------------
7865    Rmax=INT(5*C)+N-M+Rmaxbias ! (c/2r)^2=.01  See Flammer p.17.
7870  ! IF Rmax<Rmaxmin THEN Rmax=Rmaxmin
7875    IF Rmax>Rmaxmax THEN Rmax=Rmaxmax
7880    IF Lcf$="C" THEN L=FNLamda(M,N,C,Laccu)
7885    IF Lcf$="F" THEN L=FNLamdaf(M,N,C,Laccu)
7890  ! PRINT USING """L"",5X,SD.17DE";L
7895  ! --------------------Rstart------------------------------------
7900    IF N-M<-.1 THEN PRINT "ERROR N-M<0"
7905    SELECT ABS((N-M) MOD 2)
7910    CASE <.1             ! Even.
7915      Rstart=0
7920    CASE >.9             ! Odd.
7925      Rstart=1
7930    END SELECT
7935  ! -----------------------Nmr----------------------------
7940  ! Pgamma(M,Rstart,C,Pg)
7945  ! Nmr=-Pg+L
7950    Nmr=-FNGamma(M,Rstart,C)+L
7955    FOR R=Rstart TO Rmax STEP 2
7960      IF R=Rstart OR R=N-M THEN GOTO Skip1
7965      IF (N-M=0 OR N-M=1) AND (R=Rstart+2) THEN GOTO Skip1
7970      Rd=R
7975      IF R>N-M THEN Rd=R-2
7980  !   Pgamma(M,Rd,C,Pg)
7985  !   Pbc(M,Rd,Bc)
7990  !   Nmr=-Pg+L-Bc*C4/Nmr       ! By CSUB.
7995      Nmr=-FNGamma(M,Rd,C)+L-FNBc(M,Rd)*C^4/Nmr
8000  !        _________________,__________________
8005  !       Not converge when above two values are nearly equal.
8010  ! -----------------------d/d----------------------------
8015  Skip1: SELECT R-(N-M)
8020      CASE >.1
8025        Dd(R)=+Nmr/FNNd(M,R,C)    !
8030      CASE -.1 TO .1              !
8035        Dd(R)=1
8040      CASE <-.1
8045        Dd(R)=+FNNd(M,R+2,C)/Nmr  !
8050      END SELECT
8055    NEXT R
8060  ! -------------------------d----------------------------
8065    IF (N-M-2)-1<.1 THEN GOTO Skip2
8070    FOR R=N-M-4 TO Rstart STEP -2
8075      Dd(R)=Dd(R+2)*Dd(R)
8080    NEXT R
8085  Skip2: IF Rmax-(N-M+2)<.1 THEN GOTO Skip3
8090    FOR R=N-M+4 TO Rmax STEP 2
8095      Dd(R)=Dd(R-2)*Dd(R)
8100    NEXT R
8105  ! ---------------Determin Rmax--------------------------
8110  Skip3:!FOR R=Rstart TO Rmax STEP 2
8115   ! PRINT "R,D",R,Dd(R)
8120 ! NEXT R
8125    Rm=N-M
8130    Dmin=1.E+30
8135    FOR R=N-M TO Rmax STEP 2
8140      D=ABS(Dd(R))
8145      IF D<Dmin THEN
8150        Dmin=D
8155        Rm=R
8160      END IF
8165    NEXT R
8170    IF Rm=Rmax AND Dd(Rm)>Daccu THEN
8175      BEEP
8180      PRINT "Rm=Rmax in Dbyd2. Can increase accu. by Rmax. Rmax=";Rmax
8185    END IF
8190    Rmax=Rm
8195  SUBEND
8200  ! ####################################################
8205  SUB D_c
8210  ! Check program of d.
8211  ! X,4C15
8215  ! ####################################################
8220    COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
8225  In: INPUT "M,N,C(Stop=0,0,0)",M,N,C
8226    IF M+N+C=0 THEN SUBEXIT
8230    PRINT "M,N,C",M,N,C
8235    Presphe(M,N,C,L)
8240    Dminus(M,N,C,L)
8245    Kappa(M,N,C,3,K1,K2)
8250    PRINT "Kappa1,2",K1,K2
8255    FOR R=Rmin TO Rmax STEP 2
8260      PRINT USING """R="",SDDD,5X,""D="",S.4DE";R,Dd(R)
8265    NEXT R
8270    GOTO In
8275  SUBEND
8280  ! ####################################################
8285  SUB Lgna2_der(Mmax,Nmax,Nmin,X,Q(*),Qd(*))
8290  !  Derivative of associate Legendre func. of 2nd kind in matrix form.
8295  !  See Lgna_der which is not in matrix form.
8300  ! ####################################################
8305    REDIM Q(0:Mmax,Nmin:Nmax+1),Qd(0:Mmax,Nmin:Nmax)
8310    Lgna2(Mmax,Nmax+1,Nmin,X,Q(*))
8315    FOR M=0 TO Mmax
8320      FOR N=Nmin TO Nmax
8325        Qd(M,N)=((M-N-1)*Q(M,N+1)+(N+1)*X*Q(M,N))/(1-X*X)
8330      NEXT N
8335    NEXT M
8340  SUBEND
8345  SUB Abybtoxi
8350  ! a/b to xi of prolate spheroid.
8355  ! ####################################################
8360    INPUT "a/b Start, End, Step ?",S,E,Inc
8365    PRINT
8370    PRINT "a/b to xi of prolate spheroid"
8375    PRINT "*****************************"
8380    PRINT
8385    PRINT "a/b       xi"
8390    PRINT
8395    FOR Ab=S TO E STEP Inc
8400      IF Ab=1 THEN Ab=Ab+1.E-10
8405      PRINT Ab,(1-(1/Ab)^2)^(-1/2)
8410    NEXT Ab
8415  SUBEND
8420  ! ####################################################
8425  SUB Lamdac
8430    ! Check lamda.
8435    ! ####################################################
8440  ! COM /Lf/ Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup,L(6,12,240)
8445    COM /Lcf/Lcf$[8]
8450    DIM Al$[120]
8455    READ Laccu,Lmin,Lmax,Cerror
8460    DATA 1E-15,.1,1000,1E-10
8465    INPUT "Value(V), Graph1(G1:L,Ls,Ll-ka), Graph2(G2:L,n-ka) ?",Vg$
8470    ! ------------------------------------------------------
8475    SELECT Vg$
8480    CASE "V"
8485      INPUT "M,N,C (C=0 to Exit) ?",M,N,C
8490      IF C=0 THEN SUBEXIT
8495      Lpre(M,N)
8500      L=FNLamda(M,N,C,Laccu)
8505      PRINT USING Imglc;M,N,C,L
8510      BEEP
8515      GOTO 8485
8520  Imglc: IMAGE "M,N,C,L",2(3X,DD),3X,DDD.D,3X,S.DDDDDDDE
8525      ! ------------------------------------------------------
8530    CASE "G1"
8535      GOSUB Cin
8540      ALLOCATE Lmd(1:3,1:Icmax),Ll$(1:3)[10]
8545      Ll$(1)="Lamda"
8550      Ll$(2)="Lamda_sc"
8555      Ll$(3)="Lamda_lc"
8560      INPUT "One_m_n or Multi_m_n ?",Om$
8565      IF Om$="M" THEN
8570        INPUT "Sequential plotting (Y/N) ?",Sp$
8575        IF Sp$="Y" THEN CALL Seqplot
8580        INPUT "Mmin,Mmax,Nmax ?",M1,M2,N2
8585      ELSE
8590        Sp$="N"
8595        INPUT "M,N ?",M1,N2
8600        M2=M1
8605      END IF
8610      FOR M=M1 TO M2
8615        IF Om$="M" THEN
8620          N1=M
8625        ELSE
8630          N1=N2
8635        END IF
8640        FOR N=N1 TO N2
8645          PRINT
8650          PRINT "M ";M;"   N ";N
8655          PRINT
8660          PRINT "  No.       c         L           Ls          Ll"
8665          Lpre(M,N)
8670          C=Cmin
8675          Ic=1
8680          WHILE C<Cmax+Cerror
8685            Lmd(1,Ic)=FNLamda(M,N,C,Laccu)
8690            Lmd(2,Ic)=FNLamdasc(M,N,C)
8695            Lmd(3,Ic)=FNLamdalc(M,N,C)
8700            PRINT USING "DDD,7X,DD.DD,5X,3(SD.DDE,3X)";Ic,C,Lmd(1,Ic),Lmd(2,Ic),Lmd(3,Ic)
8705            FOR I=1 TO 3
8710              IF Lmd(I,Ic)<Lmin THEN Lmd(I,Ic)=Lmin
8715              IF Lmd(I,Ic)>Lmax THEN Lmd(I,Ic)=Lmax
8720            NEXT I
8725            GOSUB Linlog
8730          END WHILE
8735          BEEP 2000,3
8740          Al$="Lamdac G1 m="&VAL$(M)&" n="&VAL$(N)
8745          Plot(Lmd(*),Cmin,Cmax,Cinc,Linlog$,"c","Lamda"," ",(Sp$="Y"),Al$,Ll$(*))
8750        NEXT N
8755      NEXT M
8760      DEALLOCATE Lmd(*),Ll$(*)
8765      GOTO 8475
8770      ! ------------------------------------------------------
8775    CASE "G2"
8780      GOSUB Cin
8785      INPUT "Sequential plot (Y/N) ?",Sp$
8790      IF Sp$="Y" THEN
8795        Seqplot
8800      END IF
8805      Prelfile(Paraok)    !********
8810      INPUT "Mmin,Mmax,Nmax ?",Mmin,Mmax,Nmax
8815      FOR M=Mmin TO Mmax
8820        ALLOCATE Lmd(M:Nmax,1:Icmax),Ll$(M:Nmax)[10]
8825        FOR N=M TO Nmax
8830          PRINT
8835          PRINT "M ";M;"   N ";N
8840          PRINT
8845          PRINT "  No.       c         L "
8850          Ll$(N)="n="&VAL$(N)
8855          IF Lcf$="C" THEN CALL Lpre(M,N)
8860          C=Cmin
8865          Ic=1
8870          WHILE C<Cmax+Cerror
8875            IF Lcf$="C" THEN Lmd(N,Ic)=FNLamda(M,N,C,Laccu)
8880            IF Lcf$="F" THEN Lmd(N,Ic)=FNLamdaf(M,N,C,Laccu)
8885            PRINT USING "DDD,7X,DD.DD,5X,SD.15DE";Ic,C,Lmd(N,Ic)
8890            IF Lmd(N,Ic)<Lmin THEN Lmd(N,Ic)=Lmin
8895            GOSUB Linlog
8900          END WHILE
8905        NEXT N
8910        BEEP 2000,3
8915        Al$="Lamdac G2 m="&VAL$(M)&" nmax="&VAL$(Nmax)
8920        Plot(Lmd(*),Cmin,Cmax,Cinc,Linlog$,"c","Lamda"," ",(Sp$="Y"),Al$,Ll$(*))
8925        DEALLOCATE Lmd(*),Ll$(*)
8930      NEXT M
8935    END SELECT
8940    SUBEXIT
8945  Cin:  ! ------------------------------------------------------
8950    INPUT "Cmin,Cmax,Cinc,LIN_or_LOG ?",Cmin,Cmax,Cinc,Linlog$
8955    IF Linlog$="LOG" THEN
8960      Icmax=INT(LGT(Cmax/Cmin)/LGT(Cinc)+Cerror)+1
8965    ELSE
8970      Icmax=INT((Cmax-Cmin)/Cinc+Cerror)+1
8975    END IF
8980    RETURN
8985  Linlog: !-----------------------------------------------------
8990    IF Linlog$="LOG" THEN
8995      C=C*Cinc
9000    ELSE
9005      C=C+Cinc
9010    END IF
9015    Ic=Ic+1
9020    RETURN
9025  SUBEND
9030  ! ####################################################
9035  SUB Srmnc
9040  ! Smn,Rmn1 & Rmn2 as a functions of ka or h.
9045  ! m,n,Xi & Theta are parameters.
9046  ! X,5106
9050  ! ####################################################
9055    DEG
9060    !Prelfile(Paraok)
9065  Once: INPUT "M,N,Xi(>1),Theta ?",M,N,Xi,Theta
9066    INPUT "Is variable ka or h (K/H) ?",Kh$
9070    Eta=COS(Theta)
9071    IF Kh$="K" THEN
9075      INPUT "Ka_min, _max, _step ?",Khmin,Khmax,Khstep
9076      Xl$="ka"
9077    ELSE
9078      INPUT "h_min, _max, _step ?",Khmin,Khmax,Khstep
9079      Xl$="h"
9080    END IF
9082    PRINT
9085    PRINT "Smn, Rmn1, & Rmn2 as a function of  ka or h."
9090    PRINT "********************************************"
9095    PRINT
9100    PRINT "  m=";M;"  n=";N;"  Xi=";Xi;"  Theta(Eta)=";Theta;"(";PROUND(Eta,-2);")"
9101    IF Kh$="K" THEN
9105      PRINT "  ka_min ( _step) _max=";Kamin;" (";Kastep;")";Kamax
9106    ELSE
9107      PRINT "  h_min ( _step) _max=";Hmin;" (";Hstep;")";Hmax
9108    END IF
9110    PRINT
9115    INPUT "Parameters OK ? (Y/N)",Ok$
9120    IF Ok$="N" THEN Once
9125    Prelfile(Paraok)   !***********
9130       ! ------------------------------------------------------
9135    PRINT " ka        h        Smn         Rmn1        Rmn2         Err"
9140    Imax=INT((Khmax-Khmin)/Khstep)
9145    ALLOCATE Ans(3,Imax),Ll$(3)[10]
9150    I=0
9155    FOR Kh=Khmin TO Khmax STEP Khstep
9156      IF Kh$="K" THEN
9160        H=Kh/Xi
9161        Ka=Kh
9162      ELSE
9163        H=Kh
9164        Ka=Kh*Xi
9165      END IF
9167      Presphe(M,N,H,L)
9170      Ans(0,I)=FNSmn(M,N,H,Eta)
9175      Rmnh(3,M,N,H,L,Xi,Ans(1,I),Ans(2,I),Rmn1d,Rmn2d,Ans(3,I))
9180      PRINT USING "2(DD.DD,5X),4(SD.DDE,3X)";Ka,H,Ans(0,I),Ans(1,I),Ans(2,I),Ans(3,I)
9185      I=I+1
9190    NEXT Kh
9195    BEEP 2000,3
9200       ! ------------------Plot--------------------------------
9201    INPUT "Plot, More, Exit ? (P/M/E)",Pme$
9202    IF Pme$="M" THEN GOTO Once
9203    IF Pme$="E" THEN SUBEXIT
9205    DIM Al$[120]
9210    Al$="Srmnc m="&VAL$(M)&" n="&VAL$(N)&" Xi="&VAL$(Xi)&" Theta="&VAL$(Theta)
9215    Ll$(0)="Smn"
9220    Ll$(1)="Rmn1"
9225    Ll$(2)="Rmn2"
9230    Ll$(3)="Err"
9235  Plot: Plot(Ans(*),Khmin,Khmax,Khstep,"LIN",Xl$,"Smn,Rmn,Err"," ",0,Al$,Ll$(*),1)
9240    DEALLOCATE Ans(*),Ll$(*)
9245    GOTO Once
9250  SUBEND
9255  ! ####################################################
9260  DEF FNInte(X)
9265    !  Examine if real X is integer.
9270    !  If X is integer, return X.0, not INT(X)+1.5.
9275    ! ####################################################
9280    Xerror=1.E-10
9285    Xf=FRACT(X)
9290    Xi=INT(X)
9295    SELECT Xf
9300    CASE <Xerror
9305      RETURN Xi
9310    CASE >1-Xerror
9315      RETURN Xi+1
9320    CASE ELSE
9325      RETURN Xi+.5
9330    END SELECT
9335  FNEND
9340  ! ####################################################
9345  SUB Scatka
9350  ! Frequency characteristics of backscattering for soft, liquid & rigid spheroid.
9355  ! ####################################################
9360    DEG
9365    DIM Sf(2,0,1)     ! 0 because of one dir.
9370    PRINT "Frequency characteristics of backscattering for soft, liquid & rigid spheroid."
9375    PRINT "******************************************************************************"
9380    INPUT "Kind(0:S, 1:S,L,R, 2:S,R),Mmax,Nmax ?",Kind,Mmax,Nmax
9385    IF Kind=1 THEN
9390      INPUT "Method for liquid (Approx./Exac.) ?",Ae$
9395      IF Ae$="E" THEN ALLOCATE Sfe(0,1)
9400    END IF
9405    INPUT "Normalization ? (1:1/2,2:1/ho,3:1/koa,4:2/koa,5:2/koa*(a/b)^2,6:1/2koa)",Norm
9410    SELECT Norm
9415    CASE 1
9420      Nf$="1/2"                        ! Yeh              (kg/2)
9425    CASE 2
9430      Nf$="1/ho"                       ! Spence, Tobocman (g/q)
9435    CASE 3
9440      Nf$="1/k0a"                      ! Furusawa         (g/a)
9445    CASE 4
9450      Nf$="2/k0a"                      ! Varadan          (g/(a/2))
9455    CASE 5
9460      Nf$="2/k0a*(a/b)^2"              ! Senior, Siegel   (g/g0,squared,
9465    CASE 6                             !   Hickling          end-on)
9470      Nf$="1/2k0a"                     ! Furusawa         (g/2a,dB)
9475    END SELECT
9480    INPUT "A/B,Theta ?",Abyb,Theta
9485    Thetas=-10
9490    IF Norm=4 THEN INPUT "Thetas",Thetas
9495    INPUT "K0a: Min,Max,Step ?",K0amin,K0amax,K0astep
9500    Prelfile(Paraok,K0amin,K0amax,K0astep,Xi0,Abyb)
9505    IF Kind=1 THEN INPUT "C0byC1,Rho10 ?",C0byc1,Rho10
9510    Xi0=1/SQR(1-1/Abyb/Abyb)
9515    T1=TIMEDATE
9520  ! ----------------------F & f---------------------------
9525    I=0
9530    FOR K0a=K0amin TO K0amax STEP K0astep
9535      I=I+1
9540    NEXT K0a
9545    Imax=I-1
9550    ALLOCATE Bs(2,Imax)
9555    I=0
9560    FOR K0a=K0amin TO K0amax STEP K0astep
9565      H0=K0a/Xi0
9570      SELECT Norm
9575      CASE 1
9580        Nfp=1/2
9585      CASE 2
9590        Nfp=1/H0
9595      CASE 3
9600        Nfp=1/K0a
9605      CASE 4
9610        Nfp=2/K0a
9615      CASE 5
9620        Nfp=2/K0a*Abyb^2                  ! Senior
9625      CASE 6
9630        Nfp=1/(K0a*2)
9635      END SELECT
9640      IF Kind=1 THEN H1=H0*C0byc1  ! h=kq=2Pifq/c          I- one dir.
9645      IF Kind<>1 OR Ae$="A" THEN
9650        Scatfunslr(Kind,Mmax,Nmax,H0,H1,Rho10,Xi0,Theta,Thetas,0,180,Sf(*))
9655        ! Incident direction is (180-Theta,180)
9660      ELSE
9665        Scatfunslr(2,Mmax,Nmax,H0,H1,Rho10,Xi0,Theta,Thetas,0,180,Sf(*))
9670        Scatfunl(Mmax,Nmax,H0,H1,Rho10,Xi0,Theta,Thetas,0,180,Sfe(*))
9675        Sf(1,0,0)=Sfe(0,0)
9680      END IF
9685      FOR K=0 TO 2
9690        IF Norm=5 THEN
9695          Bs(K,I)=(Sf(K,0,0)*Nfp)^2    ! Square for comp.with Senior.
9700        ELSE
9705          Bs(K,I)=Sf(K,0,0)*Nfp
9710        END IF
9715        IF Bs(K,I)<=0 THEN Bs(K,I)=1.E-30
9720        PRINT USING """  k0a "",DD.DD,""    K "",D,""    Bs "",SD.DDE,""    T "",DD.DD";K0a,K,Bs(K,I),(TIMEDATE-T1)/3600
9725      NEXT K
9730      I=I+1
9735    NEXT K0a
9740    PRINT "Time=";PROUND((TIMEDATE-T1)/3600,-2)
9745    BEEP 2000,2
9750  ! -------------------Print------------------------------
9755    PRINT
9760    PRINT "Kind(0:S,1:L,2:R),Mmax,Nmax ",Kind,Mmax,Nmax
9765    PRINT "A/B,Theta ",Abyb,Theta
9770  ! PRINT "C0[m/s],H0max,H0step ",C0,H0max,H0step
9775    IF Kind=1 THEN PRINT "C0BYC1,Rho10 ",C0byc1,Rho10
9780    ! ------------------Plot--------------------------------
9785  P: INPUT "Plot kind (S/L/R/SR/SLR/N),  LIN or LOG ?",Pk$,Linlog$
9790    SELECT Pk$
9795    CASE "S"
9800      K=0
9805      GOSUB Replace
9810      Linlbl$(0)="SOFT"
9815    CASE "L"
9820      IF Kind=0 OR Kind=2 THEN GOTO P
9825      K=1
9830      GOSUB Replace
9835      Linlbl$(0)="LIQUID"
9840    CASE "R"
9845      IF Kind=0 THEN GOTO P
9850      K=2
9855      GOSUB Replace
9860      Linlbl$(0)="RIGID"
9865    CASE "SR"
9870      IF Kind=0 THEN GOTO P
9875      ALLOCATE Bp(1,Imax)
9880      ALLOCATE Linlbl$(1)[10]
9885      FOR K=0 TO 1
9890        FOR I=0 TO Imax
9895          Kk=K+(K=1)
9900          IF Linlog$="LOG" THEN
9905            Bp(K,I)=LGT(Bs(Kk,I))
9910            IF Norm=6 THEN Bp(K,I)=20*Bp(K,I)  ! dB
9915          ELSE
9920            Bp(K,I)=Bs(Kk,I)
9925          END IF
9930        NEXT I
9935      NEXT K
9940      Linlbl$(0)="SOFT"
9945      Linlbl$(1)="RIGID"
9950      K=20
9955    CASE "SLR"
9960      ALLOCATE Bp(2,Imax)
9965      ALLOCATE Linlbl$(2)[10]
9970      FOR K=0 TO 2
9975        FOR I=0 TO Imax
9980          IF Linlog$="LOG" THEN
9985            Bp(K,I)=LGT(Bs(K,I))
9990            IF Norm=6 THEN Bp(K,I)=20*Bp(K,I)  ! dB
9995          ELSE
10000           Bp(K,I)=Bs(K,I)
10005         END IF
10010       NEXT I
10015     NEXT K
10020     Linlbl$(0)="SOFT"
10025     Linlbl$(1)="LIQUID"
10030     Linlbl$(2)="RIGID"
10035     K=210
10040   CASE "N"
10045     SUBEXIT
10050   END SELECT
10055   ! --------------------Plot------------------------------
10060   DIM Ylbl$[30],Lilbl$(0)[10],Ttl$[50],Led$[80]
10065   Led$="K3"&" M="&VAL$(Mmax)&" N="&VAL$(Nmax)&" X="&VAL$(DROUND(Xi0,6))&" T="&VAL$(Theta)
10070   Led$=Led$&" A/B="&VAL$(Abyb)&" C01="&VAL$(C0byc1)&" R10="&VAL$(Rho10)
10075   BEEP
10080   IF Linlog$="LOG" THEN
10085     Ylbl$="log(F*"&Nf$&")"
10090     IF Norm=6 THEN Ylbl$="20log(F/2a)"
10095   ELSE
10100     Ylbl$="F*"&Nf$
10105   END IF
10110   Ttl$="k0a-norm.F"
10115   Plot(Bp(*),K0amin,K0amax,K0astep,"LIN","k0a",Ylbl$,Ttl$,0,Led$,Linlbl$(*),1)
10120   DEALLOCATE Bp(*)
10125   DEALLOCATE Linlbl$(*)
10130   GOTO P
10135  Replace:! ------------------------------------------------------
10140   ALLOCATE Bp(0,Imax)
10145   ALLOCATE Linlbl$(0)[10]
10150   FOR I=0 TO Imax
10155     IF Linlog$="LOG" THEN
10160       Bp(0,I)=LGT(Bs(K,I))
10165     ELSE
10170       Bp(0,I)=Bs(K,I)
10175     END IF
10180     IF Norm=6 THEN Bp(0,I)=20*Bp(0,I)  ! dB
10185   NEXT I
10190   RETURN
10195 SUBEND
10200 ! ####################################################
10205 SUB Scatpat
10210 ! Scattering pattern F/a or F/2a of soft,liquid(k0.NE.k1),& rigid spheroid.
10215 ! Can compare with Spence & Granger.
10216 ! X,5110
10220 ! ####################################################
10225   DIM Sf(2,1000,1),Al$[120],Rlbl$[20]   ! Sf ia redimensioned in Scatfunslr.
10230   DEG
10235   PRINT "Scattering pattern of soft, liquid, & hard spheroid."
10240   PRINT "****************************************************"
10245   ! ------------------Input & print para.-----------------
10246   INPUT "Check ? (Y/N)",Ck$
10247   IF Ck$="Y" THEN ! For pro check. Change values with "IF Ck=1 THEN"
10249     Kind=1
10250     Hk$="K"
10251     Bm$="M"
10252     K0a=12
10253     Bbya=.1
10254     Thetad=2
10255     Tplot=90
10256     C10=1.02
10257     Rho10=1.04
10258     Ok$="Y"
10259     Ck=1
10261   ELSE
10262     Ck=0
10263   END IF
10264   IF Ck=0 THEN INPUT "Kind(0:S/1:S,L,R/2:S,R),  Input h,Xi or ka,b/a (H/K) ?",Kind,Hk$
10265   PRINT "Kind(0:S/1:S,L,R/2:S,R)",Kind
10266   IF Ck=0 THEN INPUT "Bistatic(B) or Monostatic(M) ?",Bm$
10267   PRINT "Bistatic(B) or Monostatic(M)",Bm$
10270   INPUT "Mmax,Nmax ?",Mmax,Nmax
10275   PRINT "Mmax,Nmax",Mmax,Nmax
10280   Al$="Scatpat "&Bm$&",M="&VAL$(Mmax)&",N="&VAL$(Nmax)
10285   IF Hk$="H" THEN
10290     INPUT "H0,Xi0 ?",H0,Xi0
10295     K0a=H0*Xi0
10300     Bbya=SQR(1-1/Xi0^2)
10305     Al$=Al$&",Xi="&VAL$(Xi0)&",H0="&VAL$(H0)
10310   ELSE
10315     IF Ck=0 THEN INPUT "K0a, b/a ?",K0a,Bbya
10320     Xi0=1/SQR(1-Bbya^2)
10325     H0=K0a/Xi0
10330     Al$=Al$&",k0a="&VAL$(K0a)&",b/a="&VAL$(Bbya)
10335   END IF
10340   PRINT "H0,K0a",H0,K0a
10345   PRINT "b/a,Xi0",Bbya,Xi0
10350   IF Bm$="B" THEN
10355     Tplot=360
10360     INPUT "Theta0(incident),Thetad ?",Theta0,Thetad
10365     Al$=Al$&",T="&VAL$(Theta0)
10370     PRINT "Theta0(incident),Thetad",Theta0,Thetad
10375   ELSE
10380     Theta0=0
10385     IF Ck=0 THEN INPUT "Thetad, Tplot(Back:90/Back_&_foward:360) ?",Thetad,Tplot
10390     PRINT "Thetad",Thetad
10395   END IF
10400   IF Kind=1 THEN
10405     IF Ck=0 THEN INPUT "C1/0(Nearly eq. 1) ,Rho1/0 ?",C10,Rho10
10410     PRINT "C1/0,Rho1/0",C10,Rho10
10415     K1a=K0a/C10
10420     H1=K1a/Xi0
10425     PRINT "K1a,H1",K1a,H1
10430     Al$=Al$&",C10="&VAL$(C10)&",R10="&VAL$(Rho10)
10435   END IF
10440   !INPUT "Normalization: K0a(Spence) or 2K0a(dB) (1/2) ?",Nrm
10441   Nrm=2
10445   IF Ck=0 THEN INPUT "Para OK (Y/N) ?",Ok$
10450   IF Ok$="N" THEN 10264
10455   PRINT
10460   Prelfile(Paraok)
10465   PRINT
10470 ! ----------------------Cal.----------------------------
10475   Xi02=Xi0*Xi0
10480   Etap2=COS(Theta)^2
10485   IF Nrm=1 THEN
10490     Nf=K0a
10495     Rlbl$="F/a"
10500   ELSE
10505     Nf=2*K0a
10510     Rlbl$="20log(F/2a)"
10515   END IF
10520   IF Bm$="B" THEN Thetas=0
10525   IF Bm$="M" THEN Thetas=-1
10530   Scatfunslr(Kind,Mmax,Nmax,H0,H1,Rho10,Xi0,Theta0,Thetas,Thetad,180,Sf(*))
10535   MAT Sf=Sf/(Nf)
10540   ! ------------------Print results-----------------------
10541   INPUT "Printer (CRT/701)",C7$
10542   IF C7$="701" THEN PRINTER IS 701
10545   PRINT "  Normalized scat. amp. g/";Nrm;"a"
10550   PRINT
10555   PRINT "Theta","Soft(0)","Soft(180)","Liq.(0)","Liq.(180)","Hard(0)","Hard(180)"
10560   It=0
10565   FOR Theta=0 TO 180 STEP Thetad
10570     PRINT USING "3D,7X,6(SD.DDE,X)";Theta,Sf(0,It,0),Sf(0,It,1),Sf(1,It,0),Sf(1,It,1),Sf(2,It,0),Sf(2,It,1)
10575     It=It+1
10580   NEXT Theta
10585   Itmax=It-1
10586   PRINTER IS CRT
10590   BEEP 2000,3
10595 ! ------------------Plotting----------------------------
10600   IF Kind=0 THEN Kindp=0
10605   IF Kind=1 THEN Kindp=2
10610   IF Kind=2 THEN Kindp=1
10615   ALLOCATE Ll$(Kindp)[10]
10620   IF Tplot=360 THEN
10625     ALLOCATE Sf2(Kindp,2*Itmax)
10630   ELSE
10635     Itmax=INT(Itmax/2)
10640     ALLOCATE Sf2(Kindp,Itmax)
10645   END IF
10650   Ll$(0)="SOFT"
10655   Kp=0
10660   FOR K=0 TO Kindp
10665     IF K=1 AND Kind=1 THEN
10670       Ll$(K)="FLUID"
10675       Kp=1
10680     END IF
10685     IF (K=1 AND Kind=2) OR (K=2 AND Kind=1) THEN
10690       Ll$(K)="RIGID"
10695       Kp=2
10700     END IF
10705     FOR It=0 TO Itmax
10710       Sf2(K,It)=Sf(Kp,It,0)
10715     NEXT It
10720     IF Tplot=90 THEN GOTO 10740
10725     FOR It=Itmax+1 TO 2*Itmax
10730       Sf2(K,It)=Sf(Kp,2*Itmax-It,1)
10735     NEXT It
10740     IF Nrm=2 THEN
10745       FOR It=0 TO (1+(Tplot=360))*Itmax
10750         Sf2(K,It)=20*LGT(Sf2(K,It))
10755       NEXT It
10760     END IF
10765   NEXT K
10770   Plot(Sf2(*),0,Tplot,Thetad,"POL","A",Rlbl$," ",0,Al$,Ll$(*),1)    !A is dummy.
10775   DEALLOCATE Sf2(*),Ll$(*)
10780 SUBEND
10785 ! ####################################################
10790 SUB Scatfunslr(Kind,Mmax,Nmax,H0,H1,Rho1by0,Xi0,Theta0,Thetas,Thetad,Phi,Sf(*))
10795 ! Scattering function F=k0*f of soft, liquid(k0 .=. k1) & rigid spheroid.
10800 ! ------------------------------------------------------
10805 ! Kind        0        1        2
10810 !           soft     liquid   hard
10815 !           A(0)     A(1)     A(2)
10820 !  0         O                          result only for soft
10825 !  1         O        O        O        results for all
10830 !  2         O                 O        results for soft & hard
10835 ! H1,Rho10  dummy             dummy
10840 ! ------------------------------------------------------
10845 ! Theta0:Incident angle. Actual angle is 180-Theta0
10850 ! Thetas:Scattering angle for bistatic. Negative value for monostatic.
10855 ! Thetad:Increment of angle. 0 for one direction.
10860 ! Phi:Phi-Phi'. Scattering azimuthal angle measured from incident angle i.e. 180. 180 for back_ , 0 for foward_scatterring.
10865 ! Divide by k0 to obtain actual scattering function f.
10870 ! ####################################################
10875 ! COM /D/ Dd(-1000:1000),Rmin,Rmid,Rstart,Rmax  !*********
10880   ! ------------------Classify----------------------------
10885   IF Thetad<>0 THEN            ! Multi direction.
10890     IF Thetas>=0 THEN          !   Bistatic
10895       Kind$="BM"
10900       Theta1=Thetas
10905       Theta2=180
10910     ELSE                       !   Monostatic. Thetas<0.
10915       Kind$="MM"
10920       Theta1=Theta0
10925       Theta2=180-Theta0
10930     END IF
10935     Itmax=INT((Theta2-Theta1)/Thetad+.1)
10940   ELSE                         ! One direction. Thetad=0
10945     IF Thetas>=0 THEN Kind$="BO"!   Bistatic
10950     IF Thetas<0 THEN Kind$="MO"!   Monostatic
10955     Itmax=0
10960   END IF
10965 ! ----------------------DIM-----------------------------
10970   DIM A(2,1)                   ! Kind,ReIm. Expans. coeff.
10975   ALLOCATE Smn(Itmax)
10980   ALLOCATE Sm(2,Itmax,1,1)     ! Kind,Theta,Phi,ReIm. Sum by m.
10985   ALLOCATE Sn(2,Itmax,1)       ! Kind,Theta,ReIm. Sum by n.
10990   IF Kind$="MM" THEN ALLOCATE Snf(2,Itmax,1)
10995   REDIM Sf(2,Itmax,1)          ! Kind,Theta,Phi.
11000   DEG
11005   IF Mmax>Nmax THEN
11010     BEEP
11015     PRINT "Input error (Mmax>Nmax) in Scat_fun."
11020     SUBEXIT
11025   END IF
11030   ! ------------------------------------------------------
11035   MAT Sm=(0)
11040   FOR M=0 TO Mmax                                       !--M
11045     MAT Sn=(0)
11050     IF Kind$="MM" THEN MAT Snf=(0)
11055     FOR N=M TO Nmax                                     !  --N
11060       WAIT 1
11065       Presphe(M,N,H0,L)
11070   ! ------------------Smn---------------------------------
11075       DISP "Now in Scatfunslr.      M=";M,"N=";N,"H0=";H0
11080       Nmn0=FNSnorm(M,N)
11081       !PRINT "Theta0,-COS()",Theta0,-COS(Theta0)  !##########
11085       Smn0=FNSmn(M,N,H0,-COS(Theta0))  ! Actual theta=180-Theta0
11090       IF Kind$="BO" THEN Smn(0)=FNSmn(M,N,H0,COS(Thetas))
11095       IF Kind$="BM" OR Kind$="MM" THEN
11100         It=0
11105         FOR Theta=Theta1 TO Theta2 STEP Thetad
11106           !PRINT "Theta,COS()",Theta,COS(Theta)  !##########
11110           Smn(It)=FNSmn(M,N,H0,COS(Theta))
11115           It=It+1
11120         NEXT Theta
11125       END IF
11130   ! ------------------A,Hard & Soft-----------------------
11135       SELECT Kind
11140       CASE 0
11145         Rmnh(3,M,N,H0,L,Xi0,Rmn10,Rmn20,Rmn10d,Rmn20d,Err)
11150       CASE 1,2
11155         Rmnh(30,M,N,H0,L,Xi0,Rmn10,Rmn20,Rmn10d,Rmn20d,Err)
11160         Cdivid(-Rmn10d,0,Rmn10d,Rmn20d,A(2,0),A(2,1))! Hard.
11165       END SELECT
11170       Cdivid(-Rmn10,0,Rmn10,Rmn20,A(0,0),A(0,1))   ! Soft.
11175       IF Kind=0 OR Kind=2 THEN GOTO Sum
11180   ! -------------------A,Liquid----------------------------
11185       Presphe(M,N,H1,L)
11190       DISP "Now in Scatfunslr.      M=";M,"N=";N,"H1=";H1
11195       Rmnh(10,M,N,H1,L,Xi0,Rmn11,Rmn21,Rmn11d,Rmn21d,Err)
11200       C=Rho1by0*Rmn11/Rmn11d
11205       R=Rmn10-C*Rmn10d
11210       I=Rmn20-C*Rmn20d
11215       Cdivid(-R,0,R,I,A(1,0),A(1,1))
11220  Sum:! ------------------------------------------------------
11225       FOR It=0 TO Itmax                                 !    --Theta
11230         SELECT Kind$
11235         CASE "BO"
11240           Ssn=Smn0*Smn(It)/Nmn0
11245         CASE "BM"
11250           Ssn=Smn0*Smn(It)/Nmn0
11255         CASE "MO"
11260           Ssn=Smn0*Smn0/Nmn0
11265           IF (N-M) MOD 2=1 THEN Ssn=-Ssn  ! Smn(-X)=-Smn(X)
11270         CASE "MM"
11275           Ssn=Smn(It)*Smn(It)/Nmn0
11280           Ssnf=Ssn                        ! Foward
11285           IF (N-M) MOD 2=1 THEN Ssn=-Ssn
11290         END SELECT
11295         FOR Kk=0 TO 2
11300           FOR Ri=0 TO 1
11305             Sn(Kk,It,Ri)=Sn(Kk,It,Ri)+A(Kk,Ri)*Ssn
11310             IF Kind$="MM" THEN Snf(Kk,It,Ri)=Snf(Kk,It,Ri)+A(Kk,Ri)*Ssnf
11315           NEXT Ri
11320         NEXT Kk
11325       NEXT It                                           !    --Theta
11330     NEXT N                                              !  --N
11335     FOR It=0 TO Itmax
11340       FOR Ip=0 TO 1          ! Phi & Phi+180
11345         FOR Ri=0 TO 1
11350           FOR Kk=0 TO 2
11355             Sumn=Sn(Kk,It,Ri)
11360             IF Kind$="MM" AND Ip=1 THEN Sumn=Snf(Kk,It,Ri)
11365             Sm(Kk,It,Ip,Ri)=Sm(Kk,It,Ip,Ri)+Sumn*(1+(M<>0))*COS(M*(Phi+Ip*180))
11370           NEXT Kk
11375         NEXT Ri
11380       NEXT Ip
11385     NEXT It
11390   NEXT M                                                !--M
11395   ! ------------------------------------------------------
11400   FOR It=0 TO Itmax
11405     FOR Ip=0 TO 1
11410       FOR Kk=0 TO 2
11415         Sf(Kk,It,Ip)=FNCabs(Sm(Kk,It,Ip,1)*2,-Sm(Kk,It,Ip,0)*2)! *2/i
11420       NEXT Kk
11425     NEXT Ip
11430   NEXT It
11435 SUBEND
11440 ! ####################################################
11445 SUB Scatpatl
11450 ! ?/2106
11455 ! Scattering pattern of liquid spheroid.
11460 ! Can select variable parameter.
11465 !   V: Comparison with Yeh.  H0=H1.
11470 !   H: H0=H1
11475 !   B: Change minor radius b.  Comparison with Tobocman. H0<>H1.
11480 !   L: Chnge L/Lamda
11485 ! Use Scatfunl
11490 ! ####################################################
11495   DIM Sf(100,1)                ! Redimed in Scatfunl.
11500   DEG
11505   Prelfile(Paraok)
11510   PRINT
11515   PRINT "Scattering pattern of liquid spheroid."
11520   PRINT "**************************************"
11525   PRINT
11530   ! ------------------Input-------------------------------
11535   INPUT "Mmax, Nmax",Mmax,Nmax
11540   PRINT "Mmax=";Mmax;"    Nmax=";Nmax
11545   INPUT "Variable parameter (V/H/B/L) ?",Vp$
11550   PRINT "Vari. para. : ";Vp$
11555   INPUT "Rho1by0",Rho1by0
11560   PRINT "Rho1/0=";Rho1by0
11565   SELECT Vp$
11570   CASE "V"
11575     INPUT "Xi0 ?",Xi0
11580     PRINT "Xi0=";Xi0
11585     INPUT "Vn:Min,Max,Step ?",Vhbmin,Vhbmax,Vhbstep  ! Vn=h^3*ks(ks^2-1)
11590     PRINT "Vn:Min,Max,Step ",Vhbmin,Vhbmax,Vhbstep
11595     Nf$="2/k0"                 ! Yeh's norm. fact.  Result k0*F/2
11600   CASE "H"
11605     INPUT "Xi0 ?",Xi0
11610     PRINT "Xi0=";Xi0
11615     INPUT "H:Min,Max,Step ?",Vhbmin,Vhbmax,Vhbstep
11620     PRINT "H:Min,Max,Step ",Vhbmin,Vhbmax,Vhbstep
11625     Nf$="2a"                  ! Result F/2a.
11630   CASE "B"
11635     INPUT "a ?",A
11640     PRINT "a=";A
11645     INPUT "b:Min,Max,Step ?",Vhbmin,Vhbmax,Vhbstep
11650     PRINT "b:Min,Max,Step",Vhbmin,Vhbmax,Vhbstep
11655     INPUT "k0[/cm] ?",K0
11660     PRINT "k0[/cm]=";K0
11665     INPUT "c1/c0 ?",C10
11670     PRINT "c1/c0=";C10
11675     Nf$="q"                   ! Tobocman's norm. fact.  Result F/q.
11680   CASE "L"
11685     INPUT "b/a ?",Bba
11690     PRINT "b/a=";Bba
11695     Xi0=(1-Bba^2)^(-1/2)
11700     PRINT "Xi0=";Xi0
11705     INPUT "L/Lamda: Min,Max,Step ?",Vhbmin,Vhbmax,Vhbstep
11710     PRINT "L/Lamda: Min,Max,Step",Vhbmin,Vhbmax,Vhbstep
11715     INPUT "c1/c0 ?",C10
11720     PRINT "c1/c0=";C10
11725     Nf$="2a"
11730   END SELECT
11735   INPUT "Bistatic(B) or Monostatic(M) ?",Bm$
11740   PRINT "Bi_ or Mono_static: ";Bm$
11745   IF Bm$="B" THEN
11750     INPUT "Theta0(incident)",Theta0
11755     PRINT "Theta0=";Theta0
11760     Thetas=0
11765   ELSE
11770     Theta0=0
11775     Thetas=-1
11780   END IF
11785   INPUT "Thetad",Thetad
11790   PRINT "Thetad=";Thetad
11795   PRINT
11800   ! --------------------Cal.------------------------------
11805   T1=TIMEDATE
11810   Itmax=INT(180/Thetad)
11815   IF Vhbstep=0 THEN Vhbstep=1
11820   Kmax=INT((Vhbmax-Vhbmin)/Vhbstep)
11825   ALLOCATE Sf2(Kmax,2*Itmax),Ll$(Kmax)[10]
11830   FOR K=0 TO Kmax
11835     Vhb=Vhbmin+K*Vhbstep
11840     SELECT Vp$
11845     CASE "V"
11850       H0=(Vhb/Xi0/(Xi0+1)/(Xi0-1))^(1/3)
11855       H1=H0
11860       PRINT "H=";H0
11865       Ll$(K)="V="&VAL$(Vhb)
11870       Nf=2                             ! Yeh's norm. factor.
11875     CASE "H"
11880       H0=Vhb
11885       H1=H0
11890       Ll$(K)="H="&VAL$(Vhb)
11895       Nf=2*H0*Xi0                         ! =2*k0a
11900     CASE "B"
11905       Xi0=1/SQR(1-(Vhb/A)^2)
11910       PRINT "Xi0=";Xi0
11915       H0=K0*A/Xi0
11920       H1=H0/C10
11925       PRINT "H0=";H0;"   H1=";H1
11930       Ll$(K)="b="&VAL$(Vhb)
11935       Nf=H0
11940     CASE "L"
11945       K0a=Vhb*PI
11950       H0=K0a/Xi0
11955       H1=H0/C10
11960       PRINT "H0=";H0;"   H1=";H1
11965       Ll$(K)="L/Lamda="&VAL$(Vhb)
11970       Nf=2*K0a
11975     END SELECT
11980     Scatfunl(Mmax,Nmax,H0,H1,Rho1by0,Xi0,Theta0,Thetas,Thetad,180,Sf(*))
11985     ! ------------------Results-----------------------------
11990     PRINT
11995     MAT Sf=Sf/(Nf)
12000     PRINT "Theta","Sf(Phi=0)","Sf(Phi=180)"
12005     It=0
12010     FOR Theta=0 TO 180 STEP Thetad
12015       PRINT USING "DDD,7X,2(DDDD.DDDD,5X)";Theta,Sf(It,0)/2,Sf(It,1)/2
12020       It=It+1
12025     NEXT Theta
12030     T2=TIMEDATE
12035     PRINT "Time",INT((T2-T1)/60);" min"
12040     FOR It=0 TO Itmax
12045       Sf2(K,It)=Sf(It,0)
12050       IF Vp$="H" OR Vp$="L" THEN Sf2(K,It)=20*LGT(Sf2(K,It))
12055     NEXT It
12060     FOR It=Itmax+1 TO 2*Itmax
12065       Sf2(K,It)=Sf(2*Itmax-It,1)
12070       IF Vp$="H" OR Vp$="L" THEN Sf2(K,It)=20*LGT(Sf2(K,It))
12075     NEXT It
12080   NEXT K
12085   BEEP 2000,3
12090   ! ------------------Plot--------------------------------
12095   DIM Al$[120],Rlbl$[20]
12100   Al$="Scatpatl "&Bm$&" "&Vp$&" M="&VAL$(Mmax)&" N="&VAL$(Nmax)&" Vhb="&VAL$(Vhbmin)&"("&VAL$(Vhbstep)&")"&VAL$(Vhbmax)
12105   Al$=Al$&" T0="&VAL$(Theta0)&" R10="&VAL$(Rho1by0)
12110   SELECT Vp$
12115   CASE "B"
12120     Al$=Al$&" a="&VAL$(A)&" k0="&VAL$(K0)&" c10="&VAL$(C10)
12125   CASE "F","L"
12130     Al$=Al$&" Xi0="&VAL$(Xi0)&" b/a="&VAL$(Bba)&" c10="&VAL$(C10)
12135   CASE ELSE
12140     Al$=Al$&" Xi0="&VAL$(Xi0)
12145   END SELECT
12150   Rlbl$="F/"&Nf$
12155   IF Vp$="H" OR Vp$="L" THEN Rlbl$="20log("&Rlbl$&")"
12160   Plot(Sf2(*),0,180,Thetad,"POL","A",Rlbl$,"Scat. pat. of liq. spheroid",0,Al$,Ll$(*),1)  !A is dummy. 360 for forward scat.########
12165 SUBEND
12170 ! ####################################################
12175 SUB Scatl1d
12180 ! Liquid. One direction.
12185 ! ####################################################
12190   DIM Sf(0,1)
12195  In: INPUT "Bistatic(B) or Monostatic(M)",Bm$
12200   INPUT "Mmax,Nmax",Mmax,Nmax
12205   INPUT "H0,H1,XI0,Rho1by0",H0,H1,Xi0,Rho1by0
12210   INPUT "Theta0",Theta0
12215   IF Bm$="B" THEN
12220     INPUT "Thetas",Thetas
12225   ELSE
12230     Thetas=-1
12235   END IF
12240   Thetad=0
12245   Scatfunl(Mmax,Nmax,H0,H1,Rho1by0,Xi0,Theta0,Thetas,Thetad,180,Sf(*))
12250   BEEP
12255   PRINT
12260   PRINT "Bistatic(B) or Monostatic(M)",Bm$
12265   PRINT "Mmax,Nmax",Mmax,Nmax
12270   PRINT "H0,H1,XI0,Rho1by0",H0,H1,Xi0,Rho1by0
12275   PRINT "Theta0,Thetas,Thetad",Theta0,Thetas,Thetad
12280   PRINT "F",Sf(*)
12285   GOTO In
12290 SUBEND
12295 ! ####################################################
12300 SUB Scatfunl(Mmax,Nmax,H0,H1,Rho1by0,Xi0,Theta0,Thetas,Thetad,Phi,Sf(*))
12305   ! Scattering function of liquid spheroid. Yeh.
12310   ! This scat.fun. is ABS(F) of Yeh(1967).
12315   ! Theta0:Incident angle. Actual angle is 180-Theta0
12320   ! Thetas:Scattering angle for bistatic. Negative value for monostatic.
12325   ! Thetad:Increment of angle. 0 for one direction.
12330   ! Phi:Scattering azimuthal angle measured from incident angle i.e. 180.
12335   ! Divide by k0 to get actual scat.fun. g or f.
12340   ! ####################################################
12345   IF Thetad<>0 THEN            ! Multi direction.
12350     IF Thetas>=0 THEN          !   Bistatic
12355       Kind$="BM"
12360       Theta1=Thetas
12365       Theta2=180
12370     ELSE                       !   Monostatic
12375       Kind$="MM"
12380       Theta1=Theta0
12385       Theta2=180-Theta0
12390     END IF
12395     It=0
12400     FOR Theta=Theta1 TO Theta2 STEP Thetad
12405       It=It+1
12410     NEXT Theta
12415     Itmax=It-1
12420   ELSE                         ! One direction. Thetad=0
12425     IF Thetas>=0 THEN Kind$="BO"!   Bistatic
12430     IF Thetas<0 THEN Kind$="MO"!   Monostatic
12435     Itmax=0
12440   END IF
12445 ! ------------------------------------------------------
12450   ALLOCATE Smn0(Nmax),Nmn0(Nmax)
12455   ALLOCATE Smn(Nmax,Itmax)     ! N,Theta
12460   ALLOCATE Sm(Itmax,1,1)       ! Theta,Phi,ReIm. Sum by m.
12465   ALLOCATE Sn(Itmax,1)         ! Theta,ReIm. Sum by n.
12470   IF Kind$="MM" THEN
12475     ALLOCATE A(Nmax,Itmax,1)   ! N,Theta,ReIm
12480     ALLOCATE Snf(Itmax,1)      ! Theta,ReIm. Foward scat.
12485   ELSE
12490     ALLOCATE A(Nmax,0,1)       ! N,ReIm
12495   END IF
12500   REDIM Sf(Itmax,1)            ! Theta,Phi.
12505   DEG
12510   IF Mmax>Nmax THEN
12515     BEEP
12520     PRINT "Input error (Mmax>Nmax) in Scat_fun."
12525     SUBEXIT
12530   END IF
12535   ! ------------------------------------------------------
12540   MAT Sm=(0)
12545   FOR M=0 TO Mmax
12550     MAT Sn=(0)
12555     IF Kind$="MM" THEN MAT Snf=(0)
12560     Amnl(Kind$,M,Nmax,H0,H1,Theta0,Thetas,Thetad,Theta1,Theta2,Xi0,Rho1by0,A(*),Smn0(*),Smn(*),Nmn0(*))
12565     FOR N=M TO Nmax
12570       FOR It=0 TO Itmax                                 !    --Theta
12575         SELECT Kind$
12580         CASE "BO"
12585           Ssn=Smn0(N)*Smn(N,It)/Nmn0(N)
12590         CASE "BM"
12595           Ssn=Smn0(N)*Smn(N,It)/Nmn0(N)
12600         CASE "MO"
12605           Ssn=Smn0(N)*Smn0(N)/Nmn0(N)
12610           IF (N-M) MOD 2=1 THEN Ssn=-Ssn  ! Smn(-X)=-Smn(X)
12615         CASE "MM"
12620           Ssn=Smn(N,It)*Smn(N,It)/Nmn0(N)
12625           Ssnf=Ssn       ! Foward. Inc. & scat. are same direction.
12630           IF (N-M) MOD 2=1 THEN Ssn=-Ssn   ! Backward.
12635         END SELECT
12640         FOR Ri=0 TO 1
12645           IF Kind$="MM" THEN
12650             Snf(It,Ri)=Snf(It,Ri)+A(N,It,Ri)*Ssnf
12655             Sn(It,Ri)=Sn(It,Ri)+A(N,It,Ri)*Ssn
12660           ELSE
12665             Sn(It,Ri)=Sn(It,Ri)+A(N,0,Ri)*Ssn
12670           END IF
12675         NEXT Ri
12680       NEXT It
12685     NEXT N
12690     FOR It=0 TO Itmax
12695       FOR Ip=0 TO 1          ! Phi & Phi+180
12700         FOR Ri=0 TO 1
12705           Sumn=Sn(It,Ri)
12710           IF Kind$="MM" AND Ip=1 THEN Sumn=Snf(It,Ri)
12715           Sm(It,Ip,Ri)=Sm(It,Ip,Ri)+Sumn*(1+(M<>0))*COS(M*(Phi+Ip*180))
12720         NEXT Ri
12725       NEXT Ip
12730     NEXT It
12735   NEXT M
12740   ! ------------------------------------------------------
12745   FOR It=0 TO Itmax
12750     FOR Ip=0 TO 1
12755       Sf(It,Ip)=FNCabs(Sm(It,Ip,1)*2,-Sm(It,Ip,0)*2)! *2/i
12760     NEXT Ip
12765   NEXT It
12770 SUBEND
12775 ! ####################################################
12780 SUB Amnl(Kind$,M,Nmax,H0,H1,Theta0,Thetas,Thetad,Theta1,Theta2,Xi,Rho1by0,A(*),Smn0(*),Smn(*),Nmn0(*))
12785 ! Smn is used in Scatfunl.
12790 ! ####################################################
12795   COM /D/Dd(*),Rmin,Rmid,Rstart,Rmax
12800   ALLOCATE Rmn1(Nmax),Rmn2(Nmax),Rmn1d(Nmax),Rmn2d(Nmax),Rmk1(Nmax),Rmk1d(Nmax),Rmax_n(Nmax),Rmax_k(Nmax)
12805   ALLOCATE Ddn(Nmax,200),Ddk(Nmax,200),Nmk(Nmax)
12810   DEG
12815   IF Nmax<M THEN
12820     BEEP
12825     PRINT "Input error (Nmax<M) in Amn."
12830     SUBEXIT
12835   END IF
12840 ! ---------------Rmn,Rmk,Smn,Nmn------------------------
12845   FOR N=M TO Nmax           ! For H0 & N.
12850     Presphe(M,N,H0,L)
12855     DISP "In Amnl.   M=";M;"   N=";N;"   H0=";H0
12860     Nmn0(N)=FNSnorm(M,N)
12865     Smn0(N)=FNSmn(M,N,H0,-COS(Theta0))    ! Actual theta=180-Theta0
12870     Itmax=0
12875     IF Kind$="BO" THEN Smn(N,0)=FNSmn(M,N,H0,COS(Thetas))
12880     IF Kind$="BM" OR Kind$="MM" THEN
12885       It=0
12890       FOR Theta=Theta1 TO Theta2 STEP Thetad
12895         Smn(N,It)=FNSmn(M,N,H0,COS(Theta))
12900         It=It+1
12905       NEXT Theta
12910       IF Kind$="MM" THEN Itmax=It-1  ! Amn are func. of theta.
12915     END IF
12920     ! ------------------------------------------------------
12925     Rmnh(30,M,N,H0,L,Xi,Rmn1(N),Rmn2(N),Rmn1d(N),Rmn2d(N),Err) ! Z :dummy.
12930     FOR R=Rstart TO Rmax STEP 2
12935       Ddn(N,R)=Dd(R)                ! Memorized for alpha.
12940       IF H0=H1 THEN Ddk(N,R)=Dd(R)  ! Skip K loop for saving time.
12945     NEXT R
12950     Rmax_n(N)=Rmax
12955     IF H0=H1 THEN
12960       Nmk(N)=Nmn0(N)
12965       Rmk1(N)=Rmn1(N)
12970       Rmk1d(N)=Rmn1d(N)
12975       Rmax_k(N)=Rmax
12980     END IF
12985   NEXT N
12990   IF H0=H1 THEN Skip
12995   FOR K=M TO Nmax           ! For H1 & K.  K for L.
13000     Presphe(M,K,H1,L)
13005     DISP "In Amnl.   M=";M;"   N=";K;"   H1=";H1
13010     Nmk(K)=FNSnorm(M,K)
13015     Rmnh(10,M,K,H1,L,Xi,Rmk1(K),Z,Rmk1d(K),Z,Err)! Z for dummy.
13020     FOR R=Rstart TO Rmax STEP 2
13025       Ddk(K,R)=Dd(R)
13030     NEXT R
13035     Rmax_k(K)=Rmax
13040   NEXT K
13045  Skip:! ------------------Q,J,A-------------------------------
13050   M2=2*M
13055   Mm=M       ! A's for even L are not coupled with those for odd L.
13060   Rstart=0                          ! N-M=even.
13065   GOSUB Aqj
13070   Mm=M+1     ! Mm is start value of N or L. N-M=1,3,...
13075   Rstart=1                          ! N-M=odd.
13080   GOSUB Aqj
13085   SUBEXIT
13090  Aqj:! ------------------Aqj---------------------------------
13095   DISP "Now in Aqj of Amnl. M=";M
13100   I=1
13105   FOR K=Mm TO Nmax STEP 2       ! Get Ijmax, size of Q matrix.
13110     I=I+1
13115   NEXT K
13120   Ijmax=I-1
13125   FOR It=0 TO Itmax
13130     ALLOCATE Q(1:Ijmax,1:Ijmax+1,1:2),Jj(1:Ijmax,1:Ijmax,1:2),Det(1:2)
13135     I=1
13140     FOR K=Mm TO Nmax STEP 2     ! Row. L or K & I.
13145       J=1
13150       FOR N=Mm TO Nmax STEP 2   ! Column. N & J.
13155     ! ------------------Alpha-------------------------------
13160         Rmax=MIN(Rmax_n(N),Rmax_k(K))
13165         Nume=FNFact(Rstart+M2)    ! Rstart is the same for N & K.
13170         Deno=Rstart*2+M2+1
13175         Def=1
13180         S=Nume/Deno*Ddk(K,Rstart)*Ddn(N,Rstart)
13185         FOR R=Rstart+2 TO Rmax STEP 2
13190           Nume=Nume*(R+M2)*(R-1+M2)
13195           Deno=Deno+4
13200           Def=Def*R*(R-1)
13205           S=S+Nume/Deno/Def*Ddk(K,R)*Ddn(N,R)
13210         NEXT R
13215         Alpha=S*2/Nmk(K)          ! 2 is included in Nmk.
13220   ! ------------------Q & J-------------------------------
13225         IF Kind$="MM" THEN
13230           A1=Smn(N,It)/Nmn0(N)*Alpha
13235           IF (N-M) MOD 2=1 THEN A1=-A1! Incident ang. is 180-Theta.
13240         ELSE
13245           A1=Smn0(N)/Nmn0(N)*Alpha ! Smn may be 0 when Theta=0 or 90.
13250         END IF
13255         A2=Rho1by0*Rmk1(K)/Rmk1d(K)  ! Then R1 & I1=0 and Q matrix is sing.
13260         R1=A1*(Rmn1(N)-A2*Rmn1d(N))! In these ocasions make A=0.
13265         I1=A1*(Rmn2(N)-A2*Rmn2d(N))
13270         SELECT N MOD 4
13275         CASE 0             ! I^N=1
13280           Qr=R1
13285           Qi=I1
13290           Jr=-R1
13295           Ji=0
13300         CASE 1             !     i
13305           Qr=-I1
13310           Qi=R1
13315           Jr=0
13320           Ji=-R1
13325         CASE 2             !     -1
13330           Qr=-R1
13335           Qi=-I1
13340           Jr=R1
13345           Ji=0
13350         CASE 3             !     -i
13355           Qr=I1
13360           Qi=-R1
13365           Jr=0
13370           Ji=R1
13375         CASE ELSE
13380           BEEP
13385           PRINT "ERROR IN Amn."
13390           SUBEXIT
13395         END SELECT
13400         Q(I,J,1)=Qr
13405         Q(I,J,2)=Qi
13410         Jj(I,J,1)=Jr
13415         Jj(I,J,2)=Ji
13420         IF Nmax=Mm THEN
13425           Cdivid(Jj(1,1,1),Jj(1,1,2),Q(1,1,1),Q(1,1,2),A(Nmax,It,0),A(Nmax,It,1))
13430              ! If divisor=0 because of Smn=0 then A=inf. by Cdivid.
13435         END IF
13440         J=J+1
13445       NEXT N
13450   ! ------------------D-----------------------------------
13455       FOR Ri=1 TO 2
13460         S=0
13465         FOR J=1 TO Ijmax
13470           S=S+Jj(I,J,Ri)
13475         NEXT J
13480         Q(I,Ijmax+1,Ri)=S       ! D
13485       NEXT Ri
13490       I=I+1
13495     NEXT K
13500   ! ------------------A-----------------------------------
13505     Cinv(Q(*),Ijmax,1,Det(*),Ind)
13510     GOTO 13545     !************
13515     SELECT Ind
13520     CASE 2
13525       PRINT "Ind=2 in Cinv. Det=";Det(1);" ";Det(2);"i","Return A=0."
13530     CASE 3
13535       PRINT "Ind=3 in Cinv. Q matrix is singular. Return A=0."
13540     END SELECT
13545     I=1
13550     FOR N=Mm TO Nmax STEP 2
13555       FOR Ri=0 TO 1             ! Q:OPTION BASE=1, A:=0.
13560         IF Ind=2 OR Ind=3 THEN
13565           A(N,It,Ri)=0             ! When Smn=0.
13570         ELSE
13575           A(N,It,Ri)=Q(I,Ijmax+1,Ri+1)
13580         END IF
13585       NEXT Ri
13590       I=I+1
13595     NEXT N
13600     DEALLOCATE Jj(*),Q(*),Det(*)
13605   NEXT It
13610   RETURN
13615 SUBEND
13620 ! ####################################################
13625  New_udc: SUB New_udc(Char$,Array(*))
13630 ! Replace CHR$ by UDC. Called by Symbol.
13635  ! ####################################################
13640  !  This allows up to 50 new characters to be defined, each having up
13645  !  to thirty elements (rows in the array) for definition.
13650   OPTION BASE 1
13655   COM /Udc/Old_chars$,Size(*),Chars(*)
13660   IF LEN(Old_chars$)=50 THEN
13665     PRINT "User-defined Character table full."
13670   ELSE ! (still room)
13675     Pos=LEN(Old_chars$)+1
13680     Old_chars$[Pos]=Char$
13685     Size(Pos)=SIZE(Array,1)
13690     FOR Row=1 TO Size(Pos)
13695       FOR Column=1 TO 3
13700         Chars(Pos,Row,Column)=Array(Row,Column)
13705       NEXT Column
13710     NEXT Row
13715   END IF ! (room left?)
13720 SUBEND
13725 ! ######################(10)##########################
13730 SUB Bodyblad
13735 !  Compare contributions of body & bladder.
13740 !  Data are read from file.
13745 ! ####################################################
13750   COM /Pltd/Lnum,Dnum,Y(10,1000),Xmin,Xmax,Xinc,G$[8],Xlb$[40],Ylb$[40],Ttl$[50],Auto,Autled$[200],Linlbl$(20)[40],Manuled   ! <-- Pltfile
13755   DIM Autledn$[160]
13760   INPUT "ab/a, A(by exp. in db) ?",Abbya,A
13765   PRINT "ab/a, A",Abbya,A
13770   PRINT
13775   ! ------------------Soft--------------------------------
13780   PRINT "SOFT DATA"
13785   Pltfile(0)                 !
13790   ALLOCATE Yn(1:2,1:Dnum),Linlbln$(1:2)[10]
13795   FOR D=1 TO Dnum
13800     Yn(1,D)=Y(0,D-1)-40+20*LGT(Abbya)    ! Soft data are in Y(0,*)
13805     Yn(2,D)=A
13810   NEXT D
13815   Linlbln$(1)="BLADDER"
13820   Linlbln$(2)=VAL$(A)&"dB"
13825   K=1/PI/Abbya                           ! L/Lamda=k*ab/2Pi*2/(ab/a)
13830   Autledn$="SOFT:"&Autled$
13835   Plot(Yn(*),Xmin*K,Xmax*K,Xinc*K,"LIN","L/"&CHR$(151),"A[dB]"," ",0,"",Linlbln$(*),0)
13840   DEALLOCATE Yn(*),Linlbln$(*)
13845   ! ------------------Liquid------------------------------
13850   PRINT
13855   PRINT "LIQUID DATA"
13860   Pltfile(0)
13865   ALLOCATE Yn(1:1,1:Dnum),Linlbln$(1:1)[10]
13870   FOR D=1 TO Dnum
13875     Yn(1,D)=Y(1,D-1)-40                  ! Liquid data are in Yn(1,*).
13880   NEXT D
13885   Linlbln$(1)="BODY"
13890   Autledn$=Autledn$&" Liquid:"&Autled$[LEN(Autled$)-10]
13895   Plot(Yn(*),Xmin/PI,Xmax/PI,Xinc/PI,"LIN"," "," ","SPS",0,Autledn$,Linlbln$(*),1)                    ! ka=2Pi/Lamda*a=Pi*L/Lamda
13900 SUBEND
13905 ! ####################################################
13910 SUB Arrow_y(X,Y,Xe,Ye,Alength,Angle,P)
13915 ! Alength>0: Arrow,  <0: Y
13920 ! ####################################################
13925   PEN P
13930   MOVE X,Y
13935   DRAW Xe,Ye
13940   T=FNArctan(X,Y,Xe,Ye)
13945                                ! Arrow
13950   IF Alength>0 THEN CALL Wedge(Xe,Ye,P,180+T,Angle,Alength)
13955   PRINT "180+T",180+T  !**********
13960                                ! Y
13965   IF Alength<0 THEN CALL Wedge(Xe,Ye,P,T,Angle,ABS(Alength))
13970 SUBEND
13975 ! ####################################################
13980 DEF FNArctan(X1,Y1,X2,Y2)          ! -180 TO 180
13985 ! ####################################################
13990   IF X1=X2 THEN
13995     IF Y2>Y1 THEN RETURN 90
14000     IF Y2<Y1 THEN RETURN -90
14005   END IF
14010   At=ATN((Y2-Y1)/(X2-X1))
14015   IF X2>X1 THEN RETURN At
14020   IF X2<X1 THEN
14025     IF Y2>=Y1 THEN RETURN 180+At
14030     IF Y2<Y1 THEN RETURN -180+At
14035   END IF
14040 FNEND
14045 ! ####################################################
14050 SUB Wedge(X,Y,P,D,Angle,Length)          ! Original direction is +X.
14055 ! ####################################################
14060   PEN P                                  ! D rotate over PIVOT.
14065   MOVE X,Y                               ! *< Angle
14070   L=Length/COS(Angle/2)
14075 ! T=Angle/2-D+180
14080   T=D-Angle/2
14085   IDRAW +L*COS(T),L*SIN(T)
14090   PRINT "T,LcosT,LsinT",T,L*COS(T),L*SIN(T)  !********
14095 ! T=Angle/2+D-180
14100   T=D+Angle/2
14105   MOVE X,Y
14110   IDRAW +L*COS(T),+L*SIN(T)
14115   PRINT "T,LcosT,LsinT",T,L*COS(T),L*SIN(T)  !********
14120 SUBEND
14125 DEF FNBouwkamp(M,N,C,Linit,Accuracy)
14130 !  Bouwkamp's method to increase accu. of eigen value.
14135 !  See Flammer eq.(3.1.29), but it missed - sign.
14140 !  Linit is added Dlamda until Accuracy is reached or iteration
14145 !     exceeds Countmax.
14150 ! #####################(L)############################
14155   COM /C/C4
14160   DATA 80
14165   READ Countmax
14170   L=Linit
14175   Count=1                       ! Times of iteration of Bouwkamp's method.
14180  Again: DISP USING """Count & Current_accu in Bouwkamp."",5X,DD,5X,D.DDE";Count,Caccu
14185   ! ------------------1st term of denominator-------------
14190   SELECT N-M
14195   CASE <=3
14200     D1=1      ! D1 is the first term of denominator of eq.(3.1.29).
14205   CASE ELSE
14210     Bnb=1
14215     D1=1
14220     FOR R=N-M TO 2 STEP -2
14225  !    Pbc(M,R,Pbc)
14230       Bn=FNBc(M,R)*C4/FNNmrl(M,R,C,L)^2
14235       Bnb=Bnb*Bn
14240       D1=D1+Bnb
14245     NEXT R
14250   END SELECT
14255   ! ------------------2nd term of denominator-------------
14260   Bnb=1
14265   D2=0        ! D2 is the second term of denominator of eq.(3.1.29).
14270   IF Count<3 THEN CALL Rendfind(M,N-M+2,C,L,Rend)
14275     ! Rend is a func. of Lamda and may varies largely for small count.
14280   FOR R=N-M+2 TO Rend STEP 2
14285 !   Pbc(M,R,Pbc)
14290     Bn=FNNmru(M,R,C,L,Rend)^2/FNBc(M,R)/C4
14295     Bnb=Bnb*Bn
14300     D2=D2+Bnb
14305   NEXT R
14310   ! ------------------Dlamda or again---------------------
14315   Dlamda=-(FNNmrl(M,N-M+2,C,L)-FNNmru(M,N-M+2,C,L,Rend))/(D1+D2)
14320   Caccu=ABS(Dlamda/L)
14325   !PRINT "Count,L,DL",Count,L,Dlamda        !********
14330   IF Count>Countmax THEN
14335     BEEP
14340     PRINT "Accuracy is not satisfied in Bouwkamp."
14345     PRINT USING """  Accu="",SD.DDE,5X,""M="",DD,5X,""N="",DD,5X,""C="",SD.DDE";Caccu,M,N,C
14350     RETURN L+Dlamda
14355   END IF
14360   SELECT Caccu
14365   CASE <Accuracy
14370     RETURN L+Dlamda
14375   CASE ELSE                      ! L=L+Dlamda & again.
14380     L=L+Dlamda
14385     Count=Count+1
14390     GOTO Again
14395   END SELECT
14400   RETURN
14405 FNEND
14410 ! ####################################################
14415 SUB Seqbackpat
14420 ! Sequencial cal. of back scattering pattern. May be used NOF ave. of TS.
14425 ! Results are sequencially put into file.
14430 ! ####################################################
14435   COM /Lf/Mmaxf,Nmaxf,Cmin,Cmax,Cstep,Icmax,Laccup,L(6,12,240)  ! <-- Lamdaf
14440   DIM Al$[120]
14445   PRINT "Sequencial back scattering pattern of soft & liquid spheroid."
14450   PRINT "*************************************************************"
14455   PRINT
14460 ! ------------------Input & print para.-----------------
14465   INPUT "H0: min, max, step ?",H0min,H0max,H0step
14470   IF Hstep<>0 THEN
14475     Hnum=INT((H0max-H0min)/H0step+.1)+1
14480   ELSE
14485     Hnum=1
14490   END IF
14495   PRINT "H0: min, max, step, num",H0min,H0max,H0step,Hnum
14500   INPUT "b/a ?",Ba
14505   Xi0=1/SQR(1-Ba^2)
14510   PRINT "b/a, Xi0",Ba,Xi0
14515   INPUT "Rho1/0, C1/0 ?",R10,C10
14520   PRINT "Rho1/0, C1/0",R10,C10
14525   INPUT "Theta: min, inc ?",Tmin,Td
14530   Itmax=INT((180-2*Tmin)/Td+.1)
14535   PRINT "Theta: min, inc, num",Tmin,Td,Itmax
14540   INPUT "Parameter OK (Y/N) ?",Ok$
14545   IF Ok$="N" THEN GOTO 14455
14550   PRINT
14555 ! ------------------Cal.--------------------------------
14560   Prelfile(Paraok)
14565   PRINT
14570   Seqplot
14575   ALLOCATE Sf(2,Itmax,1)                   ! SLR, Theta, Phi
14580   ALLOCATE Sf2(1,Itmax),Ll$(1)[15]
14585   Ll$(0)="SOFT"
14590   Ll$(1)="LIQUID"
14595   T=TIMEDATE
14600   FOR H0=H0min TO H0max STEP H0step
14605     H1=H0/C10
14610     K0a=H0*Xi0
14615     Nf=K0a*2
14620     Mmax=INT(Nf*Ba)+2
14625     IF Mmax>Mmaxf THEN
14630       BEEP
14635       PRINT "Mmax=";Mmax;" is larger than Mmaxf=";Mmaxf;". Use Mmaxf."
14640       Mmax=Mmaxf
14645     END IF
14650     Nmax=INT(K0a/2)+2+Mmax
14655     IF Nmax>Nmaxf THEN
14660       BEEP
14665       PRINT "Nmax=";Nmax;" is larger than Nmaxf=";Nmaxf;". Use Nmaxf."
14670       Nmax=Nmaxf
14675     END IF
14680     PRINT "H0 ";H0,"H1 ";PROUND(H1,-3),"Mmax ";Mmax,"Nmax ";Nmax,"Time ";PROUND((TIMEDATE-T)/3600,-2);"h"
14685     Scatfunslr(1,Mmax,Nmax,H0,H1,R10,Xi0,Tmin,-1,Td,180,Sf(*))
14690     MAT Sf=Sf/(Nf)
14695 ! ------------------Plot--------------------------------
14700     Al$="Seqbackpat "&"M="&VAL$(Mmax)&" N="&VAL$(Nmax)&" b/a="&VAL$(Ba)&" H0="&VAL$(H0)&" R10="&VAL$(R10)&" C10="&VAL$(C10)&" Td="&VAL$(Td)
14705     FOR K=0 TO 1
14710       FOR It=0 TO Itmax
14715         Sf2(K,It)=20*LGT(Sf(K,It,0))
14720       NEXT It
14725     NEXT K
14730     Plot(Sf2(*),Tmin,180-Tmin,Td,"POL","A","20log(F/2a)"," ",1,Al$,Ll$(*),0)
14735   NEXT H0
14740 SUBEND
14745 ! ####################################################
14750 SUB Prelfile(Paraok,OPTIONAL Kamin,Kamax,Kastep,Xi,Abyb)
14755 ! Select FNLamda or FNLamdaf. By COM to Presphe & Dbyd2.
14760 ! Store Lamda by Lfile.
14765 ! Optional parameters for ka.
14770 ! Ka's is changed  to match to filed data.
14775 ! Either of Xi or Abyb may be 0.
14780 ! ####################################################
14785   COM /Lf/Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup,L(6,12,240)  ! <-- Lamdaf
14790   COM /Lcf/Lcf$[8]   ! --> Presphe, Dbyd2
14795   INPUT "Use filed Lamda (Y/N) ?",Uf$
14800   IF Uf$="N" THEN
14805     Lcf$="C"
14810     SUBEXIT
14815   END IF
14820   IF Nmax=0 THEN CALL Lfile
14825   Lcf$="F"
14830   ! ------------------------------------------------------
14835   PRINT
14840   PRINT "Lamda file.  Mmax=";Mmax;"   Nmax=";Nmax
14845   PRINT "             Cmin=";Cmin;"   Cmax=";Cmax;"   Cstep=";Cstep
14850   IF NPAR=1 THEN                  ! C
14855     INPUT "Parameters OK (Y/N) ?",Po$
14860     IF Po$="N" THEN
14865       Paraok=0
14870     ELSE
14875       Paraok=1
14880     END IF
14885     SUBEXIT
14890   END IF
14895   ! ------------------------------------------------------
14900   IF Xi=0 THEN
14905     IF Abyb<>0 AND Abyb<>1 THEN Xi=1/SQR(1-1/Abyb^2)
14910     IF Abyb=1 THEN Xi=100
14915   END IF
14920   IF Abyb=0 THEN Abyb=Xi/SQR(Xi*Xi-1)
14925   Cminp=Kamin/Xi
14930   Cmaxp=Kamax/Xi
14935   Cstepp=Kastep/Xi
14940   PRINT "Input para.  Kamin=";Kamin;"   Kamax=";Kamax;"   Kastep=";Kastep
14945   PRINT "             Cmin=";PROUND(Cminp,-3);"   Cmax=";PROUND(Cmaxp,-3);"   Cstep=";PROUND(Cstepp,-3)
14950   PRINT "             Abyb=";Abyb;"   Xi=";Xi
14955   IF Cminp<Cmin OR Cmaxp>Cmax THEN
14960     BEEP
14965     PRINT "C is out of range."
14970     Paraok=0
14975     SUBEXIT
14980   END IF
14985   ! ------------------Match-------------------------------
14990   Nstep=FNInte(Cstepp/Cstep)
14995 ! PRINT "Nstep,19611",Nstep  !*********
15000   IF Nstep=.5 THEN         ! Input step < Filed step
15005     BEEP
15010     PRINT "Cstep is too small. Change kastep to ";Cstep*Xi;". OK ?(Y/N)",Ok$
15015     IF Ok$="N" THEN
15020       Paraok=0
15025       SUBEXIT
15030     ELSE
15035       Cstepn=Cstep
15040     END IF
15045   ELSE
15050     Cstepn=INT(Nstep)*Cstep
15055   END IF
15060   Kastep=Cstepn*Xi
15065   IF Cstepn>Cminp THEN
15070     Ncmin=FNInte((Cminp-Cmin)/Cstep)
15075     Cminp=Cmin+(INT(Ncmin)+1)*Cstep
15080   ELSE
15085     Ncmin=FNInte((Cminp-Cmin)/Cstepn)
15090     IF FRACT(Ncmin)=.5 THEN
15095       Cminp=Cmin+(INT(Ncmin)+1)*Cstepn
15100     ELSE
15105       Cminp=Cmin+Ncmin*Cstepn
15110     END IF
15115   END IF
15120   Kamin=Cminp*Xi
15125   Cmaxp=Cminp+INT((Cmaxp-Cminp)/Cstepn)*Cstepn
15130   Kamax=Cmaxp*Xi
15135   PRINT "New para.  Kamin=";PROUND(Kamin,-3);"  Kamax=";PROUND(Kamax,-3);"   Kastep=";PROUND(Kastep,-3)
15140   PRINT "           Cmin=";Cminp;"  Cmax=";Cmaxp;"   Cstep=";Cstepn
15145   PRINT
15150   Paraok=1
15155 SUBEND
15160   ! ####################################################
15165 SUB Lfile
15170 ! Make Lamda file or store it in main memory.
15175 ! Lamda's are stored in L(M,N,C).
15180 ! ####################################################
15185   COM /Lf/Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup,L(6,12,240)
15190   INPUT "Make lamda file(M) or store it in main memory(S) ?",Ms$
15195   Cerror=1.E-10
15200   IF Ms$="S" THEN GOTO Store
15205 ! ------------------------------------------------------
15210   PRINT
15215   PRINT "LAMDA FILE"
15220   PRINT
15225   INPUT "Mmax(<=6),Nmax(<=12) ?",Mmax,Nmax
15230   PRINT "Mmax, Nmax      ",Mmax,Nmax
15235   INPUT "Cmin,Cmax,Cstep ?",Cmin,Cmax,Cstep
15240   Icmax=INT((Cmax-Cmin)/Cstep+Cerror)
15245   IF Icmax>240 THEN GOTO 15235
15250   PRINT "Cmin,Cmax,Cstep,Icmax",Cmin,Cmax,Cstep,Icmax
15255   INPUT "Laccu ?",Laccup
15260   PRINT "Laccu",Laccup
15265   REDIM L(Mmax,Nmax,Icmax)
15270   MAT L=(0)
15275 ! ------------------------------------------------------
15280   Filename(Filename$)
15285   PRINT "File name",Filename$
15290   Trl=(4+10)+(7*8)+((Mmax+1)*(Nmax+1)*(Icmax+1)*8)
15295   Rl=256
15300   Rn=INT(Trl/Rl)+1
15305   PRINT "Memory needed",Rn*Rl/1000;"[kBytes]"
15310   CREATE BDAT Filename$,Rn,Rl
15315   DISP "Writing on file: ";Filename$
15320   ASSIGN @File TO Filename$
15325   OUTPUT @File;Filename$,Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup
15330   ! ------------------------------------------------------
15335   T1=TIMEDATE
15340   FOR M=0 TO Mmax
15345     FOR N=M TO Nmax
15350       Lpre(M,N)
15355       Ic=0
15360       FOR C=Cmin TO Cmax+Cerror STEP Cstep
15365         L(M,N,Ic)=FNLamda(M,N,C,Laccup)
15370         PRINT USING """ M "",D,""  N "",DD,""  Ic "",DDD,""  C "",DD.DD,""  L "",SD.15DE,""  T "",DD.DD";M,N,Ic,C,L(M,N,Ic),(TIMEDATE-T1)/3600
15375         Ic=Ic+1
15380       NEXT C
15385     NEXT N
15390   NEXT M
15395   ! ------------------------------------------------------
15400   PRINT "Comp. time",PROUND((TIMEDATE-T1)/3600,-2);"h"
15405   OUTPUT @File;L(*)
15410   BEEP 2000,3
15415   PRINT "Lamda file is made."
15420   SUBEXIT
15425  Store: ! ------------------------------------------------------
15430   Filename(Filename$)
15435   ASSIGN @File TO Filename$
15440   ENTER @File;Rfilename$
15445   IF Filename$<>Rfilename$ THEN
15450     BEEP
15455     DISP "File name is not matched. Examine disc or file name & <CONT>."
15460     PAUSE
15465     GOTO 15430
15470   END IF
15475   ENTER @File;Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup
15480   PRINT
15485   PRINT "Store Lamda in main memory."
15490   PRINT "***************************"
15495   PRINT
15500   PRINT "File name",Filename$
15505   PRINT "Mmax=";Mmax,"Nmax=";Nmax
15510   PRINT "Cmin=";Cmin,"Cmax=";Cmax,"Cstep=";Cstep,"Icmax=";Icmax
15515   PRINT "Laccu=";Laccup
15520   REDIM L(Mmax,Nmax,Icmax)
15525   ENTER @File;L(*)
15530 ! PRINT L(0,3,99),L(2,5,239)  !*********8
15535   PRINT "Lamda file is stored in main memory."
15540   SUBEXIT
15545 SUBEND
15550   ! ####################################################
15555 DEF FNLamdaf(M,N,C,Laccu)
15560   !  Lamda by filed data.
15565   !  L must be in main by Lfile.
15570   !  If Cstep is matched & Laccu>Laccup then return stored data but
15575   !    if not then interpolate and/or Bouwkamp.
15580   ! ####################################################
15585   COM /Lf/Mmax,Nmax,Cmin,Cmax,Cstep,Icmax,Laccup,L(*)
15590   COM /C/C4
15595   C4=C^4
15600   REDIM L(Mmax,Nmax,Icmax)
15605   Cerror=1.E-10
15610   ! ------------------Error-------------------------------
15615   IF M>Mmax OR N>Nmax OR C>Cmax+Cerror OR C<Cmin-Cerror THEN
15620     BEEP
15625     PRINT "Parameter out of range in FNLamdaf. Return 1."
15630     RETURN 1                       !**************************
15635   END IF
15640   SELECT Laccu
15645   CASE <Laccup
15650     DISP "Laccu is increased but consume time."
15655   CASE Laccup
15660     DISP "Laccu matched"
15665   CASE >Laccup
15670     DISP "Laccu is increased for c_step_matched_data."
15675   END SELECT
15680   ! ------------------------------------------------------
15685   Icp=(C-Cmin)/Cstep
15690   Ic=INT(Icp+Cerror)
15695   IF (FRACT(Icp)<Cerror OR FRACT(Icp)>1-Cerror) AND Laccu>=Laccup THEN RETURN L(M,N,Ic)  ! Step matched.
15700   Cp=Cmin+Ic*Cstep
15705   Linit=L(M,N,Ic)+(L(M,N,Ic+1)-L(M,N,Ic))*(C-Cp)/Cstep
15710   RETURN FNBouwkamp(M,N,C,Linit,Laccu)
15715 FNEND
15720 SUB Plot(Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,Auto,OPTIONAL Autledp$,Linlblp$(*),Manuledp)
15725 ! Plotting of 2 dimensional array.
15730 ! G$   LIN   Linear scale.   X increase with step Xinc.
15735 !      SML   Semi-log scale.      ''         factor Xinc if Xinc>0.
15740 !                                            step -Xinc if Xinc<0.
15745 !      LOG   Log-log scale.       ''             ''
15750 !      POL   Polar scale.    X->T, Y->R, Xlb$ is only controler.
15755 ! Auto       0 for manual para. & 1 for automatic para by Seqplot.
15760 ! Autledp$   Automatic legend above frame.
15765 ! Linlblp$   Line label in either of for corners.
15770 !            Grouping is performed by identity of this.
15775 ! Manuledp   1 for manual legend ( 4 corners or any position ).
15780 ! Xlb$       If " " then labels are input by key.
15785 ! To superpose control by Ttl$="SPX", X are:
15790 !            S  Same scaling. Xmin,max,inc may be changed.
15795 !               Xlb$ & Ylb$=" "
15800 !            Y  New Y scaling.
15805 !            P  Plot points.
15810 ! Y(Kind,X_index), Linlbl$(Kind)
15815 ! ALLOCATE Y(K1:K2,I1:I2) & Linlbl$(K1:K2) must be declared in main.
15820 ! ####################################################
15825   COM /Pltr/Plot$[8],Sp$[8],Lastpen   ! Plot$<--Graphinit, Sp$,Lastpen-->
15830   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg  ! Upg<--Axlas,Polax, Max<--Graphinit
15835   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl
15840   COM /Lgnd/Le(1:5,1:20,1:4),Legen$(1:5,1:20)[20],Lnn(1:5),Fn,Fnmax  ! <-- Any_posi_led
15845   COM /Auto/Pltno,Tg$[8],Pen$[8],Xgs,Xts,Rsmin,Rsmax,Ygs,Yts,Axl,Line_pen$[8],Posi$[8],Mf$[8],Filename$[10]   ! <--Seqplot
15850   DIM Autled$[160],Groupe(1:10)
15851   SEPARATE ALPHA FROM GRAPHICS
15855   DEG
15860   READ Ymin,Ymax,Torg
15865  Ym: DATA 1E30,-1E+30,0
15870  Lt: DATA 1,0,4,1,6,5,5,3,8,4,4,2,7,8,4,4,6,5,5,2,8,5,7,5,4,3,6,3,5,5       ! For line type.
15875   ! ------------------Size--------------------------------
15880   IF RANK(Y)<>2 THEN
15885     BEEP
15890     PRINT "Rank of array must be 2."
15895     SUBEXIT
15900   END IF
15905   K1=BASE(Y,1)
15910   K2=K1+SIZE(Y,1)-1
15915   Lnum=K2-K1+1
15920   I1=BASE(Y,2)
15925   I2=I1+SIZE(Y,2)-1
15930   Dnum=I2-I1+1
15935   ! ------------------Optional----------------------------
15940   SELECT NPAR
15945   CASE 9
15950     Autled$=" "
15955     Llflg=0
15960     Manuled=0
15965   CASE 10
15970     Autled$=Autledp$
15975     Llflg=0
15980     Manuled=0
15985   CASE 11
15990     Autled$=Autledp$
15995     Llflg=1
16000     Manuled=0
16005   CASE 12
16010     Autled$=Autledp$
16015     Llflg=1
16020     Manuled=Manuledp
16025   END SELECT
16030   IF Llflg=1 THEN
16035     ALLOCATE Linlbl$(K1:K2)[20],Ll$[20]
16040     FOR K=K1 TO K2
16045       Linlbl$(K)=Linlblp$(K)
16050     NEXT K
16055     Ll$=Linlbl$(K1)
16060     Ig=1                           ! Groupe #
16065     MAT Groupe=(0)                 ! Num. of lines in a groupe.
16070     FOR K=K1 TO K2
16075       IF Ll$=Linlbl$(K) THEN
16080         Groupe(Ig)=Groupe(Ig)+1
16085       ELSE
16090         Ig=Ig+1
16095         Groupe(Ig)=Groupe(Ig)+1
16100         Ll$=Linlbl$(K)
16105       END IF
16110     NEXT K
16115     Igmax=Ig
16120   ELSE
16125     MAT Groupe=(1)
16130     Igmax=Lnum
16135   END IF
16140   ! ------------------Memorize on file--------------------
16145   IF Auto=0 THEN INPUT "Memorize on file (Y/N) ?",Mf$
16150   IF Mf$="N" THEN GOTO 16270
16155   IF Auto=0 THEN
16160     Filename(Filenamep$)
16165     Filename$=Filenamep$
16170   ELSE                             ! When automatic mode:
16175     Bfl=LEN(Filename$)-2*(Pltno<>0)! Basic file name length.CONT 7725
16180     Pltno=Pltno+1
16185     Pltno$=VAL$(Pltno)
16190     IF LEN(Pltno$)=1 THEN Pltno$="0"&Pltno$! 2 digit no. is added to
16195     Filename$=Filename$[1;Bfl]&Pltno$      !   file name.
16200   END IF
16205   Trl=(4+10)+(2*8)+(Lnum*Dnum*8)+(3*8)+(4+8)+(3*(4+20))+(4+80)+(Lnum*(8+4))+8
16210   Rl=256
16215   Rn=INT(Trl/Rl)+1
16220   PRINT "Memory needed",Rn*Rl/1000;"[kByte]"
16225   CREATE BDAT Filename$,Rn,Rl
16230   DISP "Writing on file: ";Filename$
16235   ASSIGN @File TO Filename$
16240   OUTPUT @File;Filename$,Lnum,Dnum
16245   IF Autled$[LEN(Autled$)-LEN(Filename$)+1;LEN(Filename$)]<>Filename$ THEN Autled$=Autled$&" "&Filename$
16250   OUTPUT @File;Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,Autled$,Linlbl$(*),Manuled
16255   ASSIGN @File TO *
16260   BEEP
16265 ! ----------------Axes & label--------------------------
16270   IF Ttl$="SPS" THEN GOTO 16345
16275   FOR K=K1 TO K2
16280     FOR I=I1 TO I2
16285       Ymax=MAX(Y(K,I),Ymax)
16290       Ymin=MIN(Y(K,I),Ymin)
16295     NEXT I
16300   NEXT K
16305  Plt: Yminc=Ymin
16310   Ymaxc=Ymax
16315   IF G$="POL" THEN
16320     Polax(Xmin,Xmax,Torg,Yminc,Ymaxc,1,Ylb$,Ttl$,Auto)! Yminc & Ymaxc is changed in SUB to Rsmin & Rsmax.
16325   ELSE
16330     Axlas(G$,Xmin,Xmax,Yminc,Ymaxc,1,Xlb$,Ylb$,Ttl$,Auto)
16335   END IF
16340 ! ---------------------Draw lines-----------------------
16345   IF Auto=0 THEN INPUT "Color, Line_type(<15), or all Bold lines(C/L/B) ?",Line_pen$
16350   Lt$="DL"                          ! DataLine. Used in Linepen.
16355   Pen=1
16360   Ig=1
16365   Sg=0
16370   GOSUB Skippen                     ! When superpose skip to lastpen.
16375   FOR K=K1 TO K2
16380     IF K=K1+Sg THEN                 ! When new groupe change pen.
16385       GOSUB Linepen
16390       Sg=Sg+Groupe(Ig)
16395       Ig=Ig+1
16400     END IF
16405     Yy=Y(K,I1)
16410     GOSUB Cut
16415     SELECT G$
16420     CASE "LIN"
16425       MOVE Xmin,Yy
16430     CASE "SML"
16435       MOVE LGT(Xmin),Yy
16440     CASE "LOG"
16445       MOVE LGT(Xmin),LGT(Yy)
16450     CASE "POL"
16455       R=Yy-Yminc
16460       MOVE R*COS(Xmin),R*SIN(Xmin)
16465     END SELECT
16470     X=Xmin
16475     FOR I=I1 TO I2
16480       Xx=X
16485       Yy=Y(K,I)
16490       GOSUB Cut
16495       SELECT G$
16500       CASE "SML"
16505         Xx=LGT(Xx)
16510       CASE "LOG"
16515         Xx=LGT(Xx)
16520         Yy=LGT(Yy)
16525       CASE "POL"
16530         R=Yy-Yminc
16535         Xx=R*COS(X)
16540         Yy=R*SIN(X)
16545       END SELECT
16550       DRAW Xx,Yy
16555       SELECT G$
16560       CASE "LIN"
16565         X=X+Xinc
16570       CASE "POL"
16575         X=X+Xinc
16580         IF X>Xmax THEN GOTO 16625
16585       CASE "LOG","SML"
16590         IF Xinc>0 THEN
16595           X=X*Xinc
16600         ELSE
16605           X=X-Xinc
16610         END IF
16615       END SELECT
16620     NEXT I
16625   NEXT K
16630   ! ------------------Line label--------------------------
16635   IF Llflg=0 OR Line_pen$="B" OR (Lnum=1 AND Ttl$<>"SPS") THEN GOTO 16935
16640   Lt$="LL"                          ! LabelLine. Used in Linepen.
16645   IF G$="POL" THEN
16650     Cfact=1
16655   ELSE
16660     Cfact=.8
16665   END IF
16670   Ycs=Csl*Yupg*Cfact
16675   Xcs=Csl*9/15*Xupg*Cfact
16680  Otherposi: IF Auto=0 THEN INPUT "Position of line label (UL,UR,BL,BR,NO) ?",Posi$
16685   IF Posi$="NO" THEN Pln$="N"
16690   IF Plot$="705" AND Posi$<>"NO" THEN INPUT "Plot line name(Y/N) ?",Pln$
16695   IF Posi$="BR" OR Posi$="UR" THEN
16700     Maxlen=0
16705     FOR K=K1 TO K2
16710       IF LEN(Linlbl$(K))>Maxlen THEN Maxlen=LEN(Linlbl$(K))
16715     NEXT K
16720     Len=Maxlen+7
16725   END IF
16730   SELECT Posi$
16735   CASE "UL"
16740     Xsp=Xsmin+Xsr/25
16745     Ysp=Ysmax-Ysr/15
16750   CASE "UR"
16755     Xsp=Xsmax-Xcs*Len
16760     Ysp=Ysmax-Ysr/15
16765   CASE "BL"
16770     Xsp=Xsmin+Xsr/25
16775     Ysp=Ysmin+Ycs*(Igmax+1)
16780   CASE "BR"
16785     Xsp=Xsmax-Xcs*Len
16790     Ysp=Ysmin+Ycs*(Igmax+1)
16795   CASE "NO"
16800     GOTO 16935
16805   END SELECT
16810   GOSUB Skippen
16815   CSIZE Csl*Cfact
16820   Pen=1
16825   Yl=Ysp+Ycs
16830   Ig=1
16835   Sg=0
16840   FOR K=K1 TO K2
16845     IF K=K1+Sg THEN
16850       Yl=Yl-Ycs
16855       GOSUB Linepen
16860       IF Line_pen$="C" THEN GOTO 16890! Only label.
16865       MOVE Xsp,Yl
16870       DRAW Xsp+Xcs*5,Yl
16875       IF Pln$="N" THEN GOTO 16905
16880       IF Line_pen$="L" THEN PEN 2
16885       LINE TYPE 1
16890       MOVE Xsp+Xcs*6*(Line_pen$<>"C"),Yl
16895       LORG 2
16900       Label(Linlbl$(K))
16905       LORG 2
16910       Ig=Ig+1
16915       Sg=Sg+Groupe(Ig)
16920     END IF
16925   NEXT K
16930 ! -------------------Legend-----------------------------
16935   IF Plot$="705" AND Line_pen$="C" THEN    ! To original pen2.
16940     BEEP
16945     PEN 0
16950     DISP "Change pen 2 if necessary and CONT."
16955     PAUSE
16960   END IF
16965   PEN 2
16970   LINE TYPE 1
16975   IF Autled$=" " THEN GOTO 17035
16980   Wl$="Y"
16985   IF Plot$="705" THEN INPUT "Need autolegend(Y/N) ?",Wl$
16990   IF Wl$="Y" THEN
16995     IF G$="POL" THEN Csl=Csl*1.2
17000     Autolegend(Autled$)
17005     IF G$="POL" THEN Csl=Csl/1.2
17010   END IF
17015   IF Manuled=1 AND Auto=0 THEN
17020     Fn=0
17025     Legend
17030   END IF
17035   !IF Plot$="701" THEN DUMP GRAPHICS #701
17036   IF Plot$="701" THEN
17037     MASS STORAGE IS ":,1406"
17039     LOADSUB ALL FROM "SCT/GDUMP_COLORED"
17040     Gdump_colored(CRT,701,"ROTATE",180,"OFF")
17041   END IF
17042   ! ---------------------For next work--------------------
17045   Lastpen=Lnum
17050   Pe$="E"
17055   IF Auto=0 THEN INPUT "Other plot or exit (P/E) ?",Pe$
17060   SELECT Pe$
17065   CASE "P"
17070     RESTORE Lt
17075     GCLEAR
17080     GOTO Plt
17085   CASE "E"
17090     SUBEXIT
17095   END SELECT
17100  Linepen: !-----------Linepen------------------------------
17105   SELECT Plot$
17110   CASE "705"
17115     SELECT Line_pen$
17120     CASE "C"
17125       PEN 0
17130       BEEP
17135       DISP "Change pen2 and CONT."     ! Color pen.
17140       PAUSE
17145       PEN 2
17150     CASE "L"
17155       READ L1,L2
17160       LINE TYPE L1,L2
17165       IF Lt$="LL" AND L1=6 THEN OUTPUT 705;"LT4,2"
17170     CASE "B"
17175       LINE TYPE 1
17180     END SELECT
17185   CASE "CRT","701"
17190     SELECT Line_pen$
17195     CASE "C"
17200       Pen=Pen+1+(Pen=3)        ! PEN 4 is green.
17205       PEN Pen                  ! R=2,Y=3,C=5,B=6,M=7
17210     CASE "L"
17215       READ L1,L2
17220       LINE TYPE L1,L2
17225     END SELECT
17230   END SELECT
17235   RETURN
17240  Skippen:! ------------------------------------------------------
17245   RESTORE Lt
17250   IF LEN(Ttl$)<3 THEN GOTO 17280
17255   IF Ttl$[1;2]="SP" THEN
17260     FOR I=1 TO Lastpen
17265       READ L1,L2
17270     NEXT I
17275   END IF
17280   RETURN
17285  Cut:! ------------------------------------------------------
17290   SELECT G$
17295   CASE "LIN","SML"
17300     IF Yy>Ysmax+10 THEN Yy=Ysmax+10     ! Prevent over flow.
17305     IF Yy<Ysmin-10 THEN Yy=Ysmin-10     ! Prevent under flow.
17310   CASE "POL"
17315     IF Yy>Ymaxc THEN Yy=Ymaxc
17320     IF Yy<Yminc THEN Yy=Yminc
17325   END SELECT
17330   RETURN
17335 SUBEND
17340 ! ####################################################
17345 SUB Axlas(Kind$,Xmin,Xmax,Ymin,Ymax,Axval,Xlb$,Ylb$,Ttl$,Auto)
17350 ! Axes, scale, and label for graph.
17355 ! Kind$  LIN:Linear,  SML:Semi_log,  LOG:Log_log.
17360 ! Axval 1: write axese value, 0 not.
17365 ! If Xlb$=" " then labels are input by key.
17370 ! ####################################################
17375   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg   ! Upg-->, Max<--Graphinit
17380   COM /Margin/Xlm,Xrm,Ybm,Ytm
17385   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl    ! --->
17390   COM /Pltr/Plot$[8],Sp$,Lastpen
17395   COM /Auto/Pltno,Tg$[8],Pen$[8],Xgs,Xts,Rsmin,Rsmax,Ygs,Yts,Axl,Line_pen$[8],Posi$[8],Mf$[8],Filename$[10]   ! <--Seqplot
17400   Xln=35                     ! Maximum character number.
17405   Yln=35
17410   Tln=45
17415   ALLOCATE Xl$[Xln+50],Yl$[Yln+50],Tytle$[Tln+50] ! 50 is margin for input
17420   CALL Graphinit(0,Auto)
17425   ! ------------------Input-------------------------------
17430   Xl$=Xlb$
17435   IF Plot$="705" AND Xlb$<>" " THEN
17440     INPUT "Want keyed label ?(Y/N)",Wl$
17445     IF Wl$="Y" THEN Xl$=" "
17450   END IF
17455   IF Xl$=" " THEN
17460     INPUT "X label ?",Xl$
17465     INPUT "Y label ?",Yl$
17470     INPUT "Tytle ?",Tytle$
17475   ELSE
17480     Yl$=Ylb$
17485     Tytle$=Ttl$
17490   END IF
17495   IF Auto=0 THEN
17500     INPUT "Tic(T) or grid(G) ?",Tg$
17505     INPUT "Do you want another pen ? (Y/N)",Pen$
17510   END IF
17515   ! ------------------Scaling-----------------------------
17520   SELECT Kind$
17525   CASE "LIN"
17530      ! Not auto, so you can specify plotting range.
17535     IF Auto=0 THEN
17540       PRINT "Xmin=";Xmin,"Xmax=";Xmax
17545       INPUT "X_scale_min, _Max, Grid_, and Tic_space ?",Xsmin,Xsmax,Xgs,Xts
17550     END IF
17555   CASE "SML","LOG"
17560     Xsmin=INT(LGT(ABS(Xmin+1.E-8)))
17565     Xsmax=INT(LGT(ABS(Xmax-1.E-8)))+1
17570   END SELECT
17575   Xsr=Xsmax-Xsmin
17580   SELECT Kind$
17585   CASE "LIN","SML"
17590     IF Auto=0 THEN
17595       PRINT "Ymin=";Ymin,"Ymax=";Ymax
17600       INPUT "Y_scale_min, _Max ?",Ysmin,Ysmax
17605     END IF
17610     Ysr=Ysmax-Ysmin
17615     IF Auto=0 THEN
17620       PRINT "Y_scale_min=";Ysmin,"Y_scale_max=";Ysmax,"Y_scale_range=";Ysr
17625       INPUT "Y_grid_&_label_space, Y_tic_space ?",Ygs,Yts
17630     END IF
17635   CASE "LOG"
17640     Ysmin=INT(LGT(ABS(Ymin+1.E-8)))
17645     Ysmax=INT(LGT(ABS(Ymax-1.E-8)))+1
17650     Ysr=Ysmax-Ysmin
17655   END SELECT
17660   WINDOW Xsmin,Xsmax,Ysmin,Ysmax
17665   ! --------------------Factors---------------------------
17670   Xupg=Xsr/Xgmax/(1-(Xlm+Xrm)/100)
17675   Yupg=Ysr/Ygmax/(1-(Ybm+Ytm)/100)
17680   Xlp=Ysmin-Ysr/15
17685   Typ=Ysmin-Ysr/6
17690   Csv=Ysr/14/Yupg                ! GDU
17695   Csl=Ysr/14/Yupg
17700   Cst=Ysr/12/Yupg
17705   Ticlen=Ysr/50/Yupg
17710   ! --------------------Axes and grid---------------------
17715   MOVE Xsmin,Ysmin
17720   DRAW Xsmax,Ysmin
17725   Ticlenx=Ticlen*Xupg
17730   Ticleny=Ticlen*Yupg
17735   IF Tg$="G" THEN Ticlenx=Xsr/2
17740   SELECT Kind$
17745   CASE "LIN","SML"
17750     FOR Decade=Ysmin TO Ysmax STEP Ygs
17755       FOR Units=1-(Decade=Ysmin) TO Ygs/Yts
17760         Y=Decade+Units*Yts
17765         MOVE Xsmin,Y
17770         DRAW Xsmin+Ticlenx*(1+(Units=Ygs/Yts)),Y
17775     !   IF Ygs<>Yts THEN
17780         MOVE Xsmax-Ticlenx*(1+(Units=Ygs/Yts)),Y
17785         DRAW Xsmax,Y
17790     !   END IF
17795       NEXT Units
17800     NEXT Decade
17805   CASE "LOG"
17810     FOR Decade=Ysmin TO Ysmax
17815       FOR Units=1 TO 1+8*(Decade<Ysmax)
17820         Y=Decade+LGT(Units)
17825         MOVE Xsmin,Y
17830         DRAW Xsmin+Ticlenx*(1+(Units=1)),Y
17835         MOVE Xsmax-Ticlenx*(1+(Units=1)),Y
17840         DRAW Xsmax,Y
17845       NEXT Units
17850     NEXT Decade
17855   END SELECT
17860   IF Tg$="G" THEN Ticleny=Ysr/2
17865   SELECT Kind$
17870   CASE "LIN"
17875     FOR Decade=Xsmin TO Xsmax STEP Xgs
17880       FOR Units=1 TO Xgs/Xts
17885         X=Decade+Units*Xts
17890         MOVE X,Ysmin
17895         DRAW X,Ysmin+Ticleny*(1+(Units=Xgs/Xts))
17900         MOVE X,Ysmax-Ticleny*(1+(Units=Xgs/Xts))
17905         DRAW X,Ysmax
17910       NEXT Units
17915     NEXT Decade
17920   CASE "SML","LOG"
17925     FOR Decade=Xsmin TO Xsmax
17930       FOR Units=1 TO 1+8*(Decade<Xsmax)
17935         X=Decade+LGT(Units)
17940         MOVE X,Ysmin
17945         DRAW X,Ysmin+Ticleny*(1+(Units=1))
17950         MOVE X,Ysmax-Ticleny*(1+(Units=1))
17955         DRAW X,Ysmax
17960       NEXT Units
17965     NEXT Decade
17970   END SELECT
17975   ! --------------------X axis value----------------------
17980   IF Pen$="Y" THEN PEN 2
17985   CLIP OFF
17990   LORG 6
17995   IF Axval=0 THEN GOTO Xl
18000   CSIZE Csv
18005   SELECT Kind$
18010   CASE "LIN"
18015     FOR I=Xsmin TO Xsmax STEP Xgs
18020       MOVE I,Ysmin
18025       LABEL USING "#,K";I
18030     NEXT I
18035   CASE "SML","LOG"
18040     FOR I=Xsmin TO Xsmax
18045       MOVE I,Ysmin
18050       LABEL USING "#,K";10^I
18055     NEXT I
18060   END SELECT
18065  Xl:! --------------------X label---------------------------
18070   CALL Label_length(Xl$,Xln)
18075   CSIZE Csl
18080   MOVE Xsmax-Xsr/2-.3*Csl*Xupg*LEN(Xl$),Xlp
18085   Label(Xl$)
18090   ! --------------------Y axis value----------------------
18095   Maxlen=1
18100   LORG 8
18105   CSIZE Csv
18110   IF Kind$="LOG" THEN Ygs=1
18115   FOR I=Ysmin TO Ysmax STEP Ygs
18120     Yav=I
18125     IF ABS(Yav)<Ygs*10^(-10) THEN Yav=0
18130     IF Kind$="LOG" THEN Yav=10^I
18135     Length=LEN(VAL$(Yav))
18140     IF Length>Maxlen THEN Maxlen=Length
18145     MOVE Xsmin,I
18150     IF Axval=1 THEN LABEL USING "#,K,X";Yav
18155   NEXT I
18160   ! --------------------Y label---------------------------
18165   CSIZE Csl
18170   LORG 4
18175   Ylpx=Xsmin-(Maxlen+1)*Csv*Xupg*.6
18180   DEG
18185   LDIR 90
18190   CALL Label_length(Yl$,Yln)
18195   MOVE Ylpx,Ysmax-Ysr/2-.3*Csl*Yupg*LEN(Yl$)
18200   Label(Yl$)
18205   ! --------------------Tytle-----------------------------
18210   LORG 6
18215   LDIR 0
18220   CSIZE Cst
18225   CALL Label_length(Tytle$,Tln)       ! Check label length.
18230   Txp=(Xsmax+Ylpx)/2-.3*Cst*Xupg*LEN(Tytle$)
18235   MOVE Txp,Typ
18240   Label(Tytle$)
18245   CLIP ON
18250 SUBEND
18255 ! ####################################################
18260 SUB Autolegend(Al$)
18265 ! Automatic legend.
18270 ! Al$ must be given such as Al$="X="&VAL$(X)&RPT$(" ",3)&"Y="&VAL$(Y)
18275 ! Character size & line number are automatically decided.
18280 ! ####################################################
18285   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl   ! <--- Axlas,Polax
18290   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg
18295   CLIP OFF
18300  Again: Length=LEN(Al$)
18305   Nl=Xsr/(Csl*Xupg*9/15)        ! Csl:GDU, Height
18310   Factor=.7
18315   Nf=Nl/Factor
18320   SELECT Length
18325   CASE <Nl
18330     Charsize=1
18335     Line=1
18340   CASE <Nf                       ! Small character.
18345     Charsize=Factor
18350     Line=1
18355   CASE <2*Nf                     ! Small character & 2 lines.
18360     Charsize=Factor
18365     Line=2
18370   CASE ELSE                      ! Lessen characters & do.
18375     Label_length(Al$,2*Nf)
18380     GOTO Again
18385     Charsize=Factor
18390     Line=2
18395   END SELECT
18400   CSIZE Csl*Charsize
18405   PEN 2
18410   Y=Ysmax+Csl*Yupg/2
18415 ! X=Xsmin
18420   LORG 4
18425   LDIR 0
18430   SELECT Line
18435   CASE 1
18440     X=Xsmin+Xsr/2-Csl*Xupg*Length*Charsize/2*(9/15)
18445     MOVE X,Y
18450     Label(Al$)
18455   CASE 2
18460     X=Xsmin+Xsr/2-Csl*Xupg*Nf*Charsize/2*(9/15)
18465     Hl=INT(Nl/Charsize)
18470     MOVE X,Y+Csl*Yupg*Charsize*1.2
18475     Label(Al$[1,Hl])
18480     MOVE X,Y
18485     LORG 4
18490     Label(Al$[Hl+1])
18495   END SELECT
18500   CLIP ON
18505 SUBEND
18510 ! ####################################################
18515 SUB Legend
18520 ! ####################################################
18525   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg
18530   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl
18535   COM /Lgnd/Le(*),Legen$(*),Lnn(*),Fn,Fnmax   ! See Any_posi_led.
18540   !                   Fn is reset in Plot & counted in this SUB.
18545   COM /Pltr/Plot$,Sp$,Lastpen
18550   DIM L$[35]
18555   Labelpen=2
18560   IF Plot$="705" THEN
18565     GOSUB Plot705
18570     SUBEXIT
18575   END IF
18580  Otherposi: INPUT "Legend position? (UL,UR,BL,BR,ANY,NO)",Posi$
18585   IF Posi$<>"NO" THEN
18590     Fn=Fn+1
18595     Fnmax=Fn
18600   END IF
18605   SELECT Posi$
18610   CASE "UL"
18615     Xsp=Xsmin+Xsr/25
18620     Ysp=Ysmax-Ysr/15
18625   CASE "UR"
18630     Xsp=Xsmax-Csl*Xupg*8
18635     Ysp=Ysmax-Ysr/15
18640   CASE "BL"
18645     Xsp=Xsmin+Xsr/25
18650     Ysp=Ysmin+Csl*Yupg*5
18655   CASE "BR"
18660     Xsp=Xsmax-Csl*Xupg*8
18665     Ysp=Ysmin+Csl*Yupg*5
18670   CASE "ANY"
18671     CONTROL KBD,15;1    !Emulate series 200
18672     CONTROL CRT,12;0
18673     LOAD KEY
18675     CALL Any_posi_led
18676     CONTROL KBD,15;0    !Return series 300
18678     CONTROL CRT,12;2
18679     LOAD KEY
18680   CASE "NO"
18685     SUBEXIT
18690   END SELECT
18695   CSIZE Csl
18700   Y=Ysp+Csl*Yupg
18705   Ln=1
18710   PEN Labelpen
18715  Next_legend: LINPUT "Legend ?   ( No more legend: ""NO""  Any other position: ""OP"" )",L$
18720   IF L$="NO" THEN
18725     Lnn(Fn)=Ln-1
18730     SUBEXIT
18735   END IF
18740   IF L$="OP" THEN
18745     Lnn(Fn)=Ln-1
18750     GOSUB Otherposi
18755   END IF
18760   Y=Y-Csl*Yupg
18765   Le(Fn,Ln,1)=Xsp
18770   Le(Fn,Ln,2)=Y
18775   MOVE Xsp,Y
18780   LORG 2
18785   LDIR 0
18790   Label(L$)
18795   Ln=Ln+1
18800   Le(Fn,Ln,3)=2
18805   Legen$(Fn,Ln)=L$
18810   Le(Fn,Ln,4)=0
18815   GOTO Next_legend
18820  Plot705: ! ------------------------------------------------------
18825   FOR Fn=1 TO Fnmax
18830     FOR Ln=1 TO Lnn(Fn)
18835       GOSUB Drawlabel
18840     NEXT Ln
18845   NEXT Fn
18850   RETURN
18855  Drawlabel:! ------------------------------------------------------
18860   PEN Labelpen
18865   LINE TYPE 1
18870   MOVE Le(Fn,Ln,1),Le(Fn,Ln,2)
18875   SELECT Legen$(Fn,Ln)
18880   CASE "SlD"
18885     DRAW Le(Fn,Ln,3),Le(Fn,Ln,4)
18890   CASE "DoT"
18895     LINE TYPE 4,1
18900     DRAW Le(Fn,Ln,3),Le(Fn,Ln,4)
18905   CASE "ArW"
18910     Arrowaniso(Le(Fn,Ln,1),Le(Fn,Ln,2),Le(Fn,Ln,3),Le(Fn,Ln,4),Xer/30,30,Labelpen)
18915   CASE ELSE
18920     LORG Le(Fn,Ln,3)
18925     LDIR Le(Fn,Ln,4)
18930     Label(Legen$(Fn,Ln))
18935   END SELECT
18940   LINE TYPE 1
18945   RETURN
18950 SUBEND
18955 ! #######################(15)#########################
18960 SUB Seqplot
18965 ! Give parameter for sequencial plotting.
18970 ! Can use automatic parameter setting.
18975 ! Para. are passed by COM /Auto/, /Range/ & /Pltr/ to Plot, Axlas, Polax &13434 !   Graphinit.
18980 ! In these SUB's must be Auto=1.
18985 ! ####################################################
18990   COM /Pltr/Plot$[8],Sp$,Lastpen       ! Plot$<--Graphinit, Sp$,Lastpen<--
18995   COM /Auto/Pltno,Tg$,Pen$,Xgs,Xts,Rsmin,Rsmax,Ygs,Yts,Axl,Line_pen$,Posi$,Mf$,Filename$   ! -->
19000   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl
19005   Pltno=0                   ! Reset plot number, counted in Plot.
19010   Pen$="Y"
19015   INPUT "Plotter (CRT,701) ?",Plot$
19020   INPUT "Linear, Semilog, LoGlog, Polar graph (L/S/G/P) ?",Po$
19025   INPUT "Color, Line_type(<15) or all Bold lines (C/L/B) ?",Line_pen$
19030   IF Plot$<>"CRT" AND Line_pen$="C" THEN GOTO 19025
19035   IF Line_pen$<>"B" THEN INPUT "Position of line labels (UL,UR,BL,BR,NO) ?",Posi$
19040   INPUT "Tic or Grid (T/G) ?",Tg$
19045   INPUT "Memorize on file (Y/N) ? (If N then cannot change plot para.)",Mf$
19050   IF Mf$="Y" THEN
19055     CALL Filename(Filenamep$)! Filename$ is commoned.
19060     Filename$=Filenamep$
19065   END IF
19070   SELECT Po$
19075   CASE "L"
19080     INPUT "X_scale_min, _Max, Grid_ & Tic_space ? ",Xsmin,Xsmax,Xgs,Xts
19085     INPUT "Y_scale_min, _Max, Grid_ & Tic_space ? ",Ysmin,Ysmax,Ygs,Yts
19090   CASE "S"
19095     INPUT "Y_scale_min, _Max, Grid_ & Tic_space ? ",Ysmin,Ysmax,Ygs,Yts
19100   CASE "P"
19105     INPUT "Angle_grid_,  _Tic_space  & Axis with label ( NO:0,0,360) ?",Xgs,Xts,Axl
19110     INPUT "R_scale_min, _Max, Grid_ & Tic_space ? ",Rsmin,Rsmax,Ygs,Yts
19115   END SELECT
19120 SUBEND
19125   ! ####################################################
19130 SUB Graphinit(Iso,Auto)
19135 ! Select plotter & VIEPORT.
19140 ! Iso: 1=isotropic, 0=anisotropic.
19145 ! ####################################################
19150   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg    ! Upg<--, Max-->
19155   COM /Margin/Xlm,Xrm,Ybm,Ytm       ! -->
19160   COM /Pltr/Plot$,Sp$,Lastpen       ! Plot$-->, Sp$<--
19165   COM /Auto/Pltno,Tg$,Pen$,Xgs,Xts,Rsmin,Rsmax,Ygs,Yts,Axl,Line_pen$,Posi$,Mf$,Filename$   ! <--
19166   Fpen=11
19170   BEEP
19175   DISP "Caution! X switch OK ?"
19180   WAIT .5
19185   IF Auto=0 THEN INPUT "Plotter(CRT/705/701) ?",Plot$
19190   IF Sp$="Y" THEN
19195     PEN Fpen
19200     GOTO Margin
19205   END IF
19210   GINIT
19215   SELECT Plot$
19220   CASE "CRT","701"
19225     PLOTTER IS CRT,"INTERNAL";COLOR MAP  !###########
19230   CASE "705"
19235     PLOTTER IS 705,"HPGL"
19240     OUTPUT 705;"VS10"
19245   END SELECT
19250   GRAPHICS ON
19255   ! ------------------------------------------------------
19260   Xgmax=100*MAX(1,RATIO)
19265   Ygmax=100*MAX(1,1/RATIO)
19270  Margin: PRINTER IS CRT
19275   IF Iso=0 THEN
19280     Xlm=20
19285     Xrm=10
19290     Ybm=25
19295     Ytm=10
19300   ELSE
19305     Xlm=5
19310     Xrm=5
19315     Ybm=11
19320     Ytm=8
19325     Xby=(100-Xlm-Xrm)/(100-Ybm-Ytm)
19330   END IF
19335   IF Auto=0 THEN
19340     PRINT "Margines X_left, right, Y_bottom, top in %.  "
19345     PRINT "        Default:",Xlm,Xrm,Ybm,Ytm
19350     PRINT "        Paper(Aniso):       35,        35,        35,        35"
19355     PRINT "        OHP(Aniso):         25,        25,        25,        25"
19360   END IF
19365  Frame: VIEWPORT Xlm/100*Xgmax,(100-Xrm)/100*Xgmax,Ybm/100*Ygmax,(100-Ytm)/100*Ygmax
19370   IF Sp$="N" THEN GCLEAR
19375   IF Plot$="CRT" OR Plot$="701" THEN
19376     PEN Fpen
19377     FRAME
19378   END IF
19380   IF Auto=0 THEN
19385     INPUT "Frame OK ? (Y/N)",Frame$
19390     IF Frame$="N" THEN
19391       PEN -Fpen
19392       FRAME
19393       PEN Fpen
19395       IF Iso=0 THEN
19400         INPUT "X-left, X-right, Y-bottom, Y-top margin ?",Xlm,Xrm,Ybm,Ytm
19405         PRINT "Xl,Xr,Yb,Yt  ";Xlm;Xrm;Ybm;Ytm
19410       ELSE
19415         INPUT "X_left, right, Y_bottom margin ?",Xlm,Xrm,Ybm
19420         Ytm=100-Ybm-(100-Xlm-Xrm)/Xby! X_y ratio is determined by default.
19425       END IF
19430       GOTO Frame
19435     END IF
19440   END IF
19445   IF Plot$="705" AND Iso=0 THEN
19450     FOR F=1 TO 2
19455       FRAME
19460     NEXT F
19465   END IF
19470 ! INPUT "Use special symbol ?(Y/N)",Ss$
19475 ! IF Ss$="Y" THEN CALL Symbol
19480 SUBEND
19485 ! ####################################################
19490 SUB Filename(Filename$)
19495 ! Input file name, purge file.
19496 ! X,4C15
19500 ! ####################################################
19505   BEEP
19510   DISP "Insert disc and make sure MSI.   MSI OK (Y/N)",
19515   INPUT Mok$
19520   IF Mok$="N" THEN
19525     INPUT "Msi:internal_FD(FD), _HD(HD), 9122_r(21) or _l(20) ?",Un$
19530     SELECT Un$
19535     CASE "FD"
19540       MASS STORAGE IS ":,1400"
19545     CASE "HD"
19550       MASS STORAGE IS ":,1406"
19551       CAT
19553       INPUT "Directory ? (Root=R)",Dt$
19556       IF Dt$<>"R" THEN MASS STORAGE IS Dt$
19557     CASE "21"
19560       MASS STORAGE IS ":,700,1"
19565     CASE "20"
19570       MASS STORAGE IS ":,700,0"
19575     END SELECT
19580   END IF
19585   CAT
19590   INPUT "File_name(Ex:CCYMMDDC,auto<=8,else<=10) or Purge(P) or MSI(M) ?",Filename$
19595   IF Filename$="M" THEN 19510
19600   IF Filename$="P" THEN
19605     INPUT "File name (<=10) ?",Filename$
19610     DISP "Purge ";Filename$&" ",
19615     INPUT "OK (Y/N) ?   LAST !!",Pok$
19620     IF Pok$="N" THEN 19590
19625     PURGE Filename$
19630     INPUT "More purge (Y/N) ?",Mp$
19635     IF Mp$="Y" THEN GOTO 19605
19640     DISP "Delete ? char. from ";Filename$;" as new name (Ex. Auto=2,Same=0,New=10) ",
19645     INPUT Cn
19650     IF Cn=10 THEN GOTO 19590
19655     Filename$=Filename$[1;LEN(Filename$)-Cn]
19660     DISP "New file name is ";Filename$;",  OK(Y/N)",
19665     INPUT Fok$
19670     IF Fok$="N" THEN GOTO 19590
19675   END IF
19680   OFF ERROR
19685   PRINT
19690 SUBEND
19695 ! ####################################################
19700 SUB Arrowaniso(X,Y,Xa,Ya,Arlen,Arang,Pen)
19705 ! Draw arrow with length Arlen in anisotropic scale.
19710 ! ####################################################
19715   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl
19720   DEG
19725   PEN Pen
19730   Yx=Ysr/Xsr
19735   MOVE X,Y
19740   DRAW Xa,Ya
19745   T=FNArctan(X,Y/Yx,Xa,Ya/Yx)
19750   L=Arlen/COS(Arang/2)
19755   Tu=T+180-Arang/2
19760   Tl=Tu+Arang
19765   IDRAW L*COS(Tu),L*SIN(Tu)*Yx
19770   MOVE Xa,Ya
19775   IDRAW L*COS(Tl),L*SIN(Tl)*Yx
19780 SUBEND
19785   ! ####################################################
19790 SUB Any_posi_led
19795 ! Draw label, line, arrow, points at any position.
19800 ! ####################################################
19805   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl
19810   COM /Lgnd/Le(1:5,1:20,1:4),Legen$(1:5,1:20)[20],Lnn(1:5),Fn,Fnmax
19815   !          Le(Fignum, Lednum, X,Y,Org,Dir)
19820   !          Legen$: legend,    Lnn: Max.legend number,    Fn: figure #
19825   COM /Pltr/Plot$,Sp$,Lastpen
19830   DEG
19835   Symbol(0)
19840   Carsorpen=1
19845   Labelpen=2
19850   CSIZE Csl
19855   IF Plot$="705" THEN GOSUB Plot705
19860   Dir$="H"
19865   Y0=Ysmin+Ysr/2
19870   X0=Xsmin+Xsr/2
19875   Y=Y0
19880   X=X0
19885   Ln=1                               ! Legend #
19890   Flg=0                              ! Flag to indicate start/end point.
19895   ON KNOB .15 GOSUB Knob_turned
19900   ON KEY 0 LABEL "UP-DOWN" GOSUB Up_down
19905   ON KEY 5 LABEL "LEFT-RIGHT" GOSUB Left_right
19910   ON KEY 6 LABEL "EXIT" RECOVER Exit
19915   OUTPUT 1;"Find position by soft key and knob. You can exit by ""EXIT""."
19920   GOSUB Plot_carsor
19925  Idle: GOTO Idle
19930   STOP
19935  Up_down:!
19940   Dir$="V"
19945   RETURN
19950  Left_right:!
19955   Dir$="H"
19960   RETURN
19965  Knob_turned:!
19970   GOSUB Erase_carsor
19975   Count=KNOBX
19980   IF Dir$="V" THEN Y=Y+Count*Ysr/300
19985   IF Y<Ysmin THEN Y=Ysmin
19990   IF Y>Ysmax THEN Y=Ysmax
19995   IF Dir$="H" THEN X=X+Count*Xsr/300
20000   IF X<Xsmin THEN X=Xsmin
20005   IF X>Xsmax THEN X=Xsmax
20010   DISP "X ";X,"Y ";Y
20015   GOSUB Plot_carsor
20020   ON KEY 1 LABEL "POSITION OK" GOSUB Position_ok
20025   RETURN
20030  Position_ok: !
20035   IF Flg=0 THEN INPUT "Label, Figure (L/F) ?",Lf$
20040   SELECT Lf$
20045   CASE "L"
20050     INPUT "Label, Direction ?",Legen$(Fn,Ln),Ldir
20055     Lorg=5
20060     Le(Fn,Ln,1)=X
20065     Le(Fn,Ln,2)=Y
20070     Le(Fn,Ln,3)=Lorg
20075     Le(Fn,Ln,4)=Ldir
20080     GOSUB Erase_carsor
20085     GOSUB Drawlabel
20090     LDIR 0
20095   CASE "F"
20100     IF Flg=0 THEN              ! Starting point.
20105       Le(Fn,Ln,1)=X
20110       Le(Fn,Ln,2)=Y
20115       INPUT "Figure name (SlD/DoT/ChN/ArW) ?",Legen$(Fn,Ln)
20120       DISP "Decide end point by knob."
20125       Flg=1                          ! Next for end point.
20130     ELSE                             ! End point.
20135       Le(Fn,Ln,3)=X
20140       Le(Fn,Ln,4)=Y
20145       GOSUB Erase_carsor
20150       GOSUB Drawlabel
20155       Flg=0
20160     END IF
20165   END SELECT
20170   Lnn(Fn)=Ln
20175   IF Flg=0 THEN Ln=Ln+1
20180   RETURN
20185  Plot_carsor: !
20190   LORG 5
20195   PEN Carsorpen
20200   MOVE X,Y
20205   LABEL "X"
20210   RETURN
20215  Erase_carsor:!
20220   LORG 5
20225   PEN -Carsorpen
20230   MOVE X,Y
20235   LABEL "X"
20240   RETURN
20245  Drawlabel:!
20250   PEN Labelpen
20255   LINE TYPE 1
20260   MOVE Le(Fn,Ln,1),Le(Fn,Ln,2)
20265   SELECT Legen$(Fn,Ln)
20270   CASE "SlD"
20275     DRAW Le(Fn,Ln,3),Le(Fn,Ln,4)
20280   CASE "DoT"
20285     LINE TYPE 4,1
20290     DRAW Le(Fn,Ln,3),Le(Fn,Ln,4)
20295   CASE "ChN"
20300     LINE TYPE 6,5
20305     DRAW Le(Fn,Ln,3),Le(Fn,Ln,4)
20310   CASE "ArW"
20315     Arrowaniso(Le(Fn,Ln,1),Le(Fn,Ln,2),Le(Fn,Ln,3),Le(Fn,Ln,4),Xsr/30,30,Labelpen)
20320   CASE ELSE
20325     LORG Le(Fn,Ln,3)
20330     LDIR Le(Fn,Ln,4)
20335     Label(Legen$(Fn,Ln))
20340   END SELECT
20345   LINE TYPE 1
20350   RETURN
20355  Plot705:!
20360   FOR Ln=1 TO Lnn(Fn)
20365     GOSUB Drawlabel
20370   NEXT Ln
20375   SUBEXIT
20380  Exit: SUBEND
20385 ! ######################(16)##########################
20390 SUB Symbol(Cm)
20395 ! Create UDC.  Cm=1 for creation, =0 for only menu.
20400  ! ####################################################
20405   IF Cm=0 THEN
20410     GOSUB Menu
20415     SUBEXIT
20420   END IF
20425   OPTION BASE 1
20430   COM /Udc/Old_chars$,Size(*),Chars(*)
20435   REAL Sigma(7,3),Infinity(16,3),Arrowb(9,3),Arrow3(5,3),Arrow8(5,3),Box(12,3)
20440   REAL Deg(9,3),Degc(19,3),Permill(33,3),Delta(5,3)
20445   REAL Solid(3,3),Dotted(11,3),Chain(11,3)
20450   REAL Psil(18,3),Phi(13,3),Tau(8,3),Theta(12,3),Lamda(6,3)
20455   REAL Omegal(15,3),Eta(10,3),Pai(10,3),Intgrl(7,3)
20460   REAL Alpha(16,3),Rho(9,3)
20465   REAL Sbs(10,3)
20470   REAL S0(9,3),S1(5,3)
20475   REAL Circle(13,3),Cross(4,3),Tangl(4,3),Square(5,3),Dot(21,3)
20480   ! ------------------Read--------------------------------
20485   READ Sigma(*),Infinity(*),Arrowb(*),Arrow3(*),Arrow8(*),Box(*)
20490   READ Alpha(*),Rho(*),Deg(*),Degc(*),Permill(*),Delta(*)
20495   READ Solid(*),Dotted(*),Chain(*)
20500   READ Psil(*),Phi(*),Tau(*),Theta(*),Lamda(*)
20505   READ Omegal(*),Eta(*),Pai(*),Intgrl(*)
20510   READ Sbs(*)
20515   READ S0(*),S1(*)
20520   READ Circle(*),Cross(*),Tangl(*),Square(*),Dot(*)
20525   ! ------------------Data--------------------------------
20530  Sigma: DATA 7,5,-2,        7,4,-1,        1,4,-1,        5.5,8.5,-1
20535   DATA 1,13,-1,       7,13,-1,       7,12,-1
20540  Infinity: DATA 4,9,-2,        5,10,-1,       6,10,-1,       7,9,-1
20545   DATA 7,8,-1,        6,7,-1,        5,7,-1,        4,8,-1
20550   DATA 4,9,-1,        3,10,-1,       2,10,-1,       1,9,-1
20555   DATA 1,8,-1,        2,7,-1,        3,7,-1,        4,8,-1
20560  Arrowb: DATA 0,0,6,         4,4,-2,        7,8,-1,        4,12,-1
20565   DATA 4,10,-1,       1,10,-1,       1,6,-1,        4,6,-1
20570   DATA 0,0,7
20575  Arrow3: DATA 0,8,-2,   27,8,-1,       18,10,-1,      18,6,-2
20580   DATA 27,8,-1
20585  Arrow8: DATA 0,8,-2,   72,8,-1,       63,10,-1,      63,6,-2,
20590   DATA 72,8,-1
20595  Box: DATA 0,0,6,         3,0,-2,        27,0,-1,       27,15,-1
20600   DATA 0,15,-1,       0,0,-1,        3,0,-1,        3,3,-1
20605   DATA 24,3,-1,       24,12,-1,      3,12,-1,       0,0,7
20610  Alpha: DATA 6,11,-2,       6,11,-1,       6,10,-1,       5,6,-1,
20615   DATA 4.5,5,-1,      3,4,-1,        2,4,-1,        1,5,-1,
20620   DATA 1,8,-1,        2,9.5,-1,      3,10,-1,       4,10,-1,
20625   DATA 5,9,-1,        5.5,6,-1,      6,5,-1,        7,4,-1
20630  Rho: DATA 1,1,-2,      1,8,-1,        2,10,-1,       5,10,-1
20635   DATA 6,8,-1,        6,6,-1,        5,4,-1,        2,4,-1
20640   DATA 1,6,-1
20645  Deg: DATA 2,11,-2,     3,12,-1,       4,12,-1,       5,11,-1
20650   DATA 5,10,-1,       4,9,-1,        3,9,-1,        2,10,-1
20655   DATA 2,11,-1
20660  Degc: DATA 3,12,-2,       3,12,-1,       2.7,11.3,-1,   2,11,-1,
20665   DATA 1.3,11.3,-1,   1,12,-1,       1.3,12.7,-1,   2,13,-1,
20670   DATA 2.7,12.7,-1,   3,12,-1,       8,11,-2,       8,11,-1,
20675   DATA 7,12,-1,       4,12,-1,       3,10,-1,       3,6,-1,
20680   DATA 4,4,-1,        7,4,-1,        8,5,-1
20685  Permil: DATA 3,11,-2,   3,11,-1,       3,10,-1,       2.5,9,-1,
20690   DATA 1.5,9,-1,      1,10,-1,       1,11,-1,       1.5,12,-1,
20695   DATA 2.5,12,-1,     3,11,-1,       6,12,-2,       6,12,-1,
20700   DATA 1,4,-1,        5,6,-2,        5,6,-1,        5,5,-1,
20705   DATA 4.5,4,-1,      3.5,4,-1,      3,5,-1,        3,6,-1,
20710   DATA 3.5,7,-1,      4.5,7,-1,      5,6,-1,        8,6,-2,
20715   DATA 8,6,-1,        8,5,-1,        7.5,4,-1,      6.5,4,-1,
20720   DATA 6,5,-1,        6,6,-1,        6.5,7,-1,      7.5,7,-1,
20725   DATA 8,6,-1
20730  Delta: DATA 1,4,-2,  1,4,-1,        5,11,-1,       7,4,-1,
20735   DATA 1,4,-1
20740  Solid: DATA 0,8,-2,     0,8,-1,        27,8,-1
20745  Dotted: DATA 0,8,-2,    0,8,-1,        3,8,-1,        6,8,-2,
20750   DATA 9,8,-1,        12,8,-2,       15,8,-1,       18,8,-2,
20755   DATA 21,8,-1,       24,8,-2,       27,8,-1
20760  Chain: DATA 0,8,-2,     0,8,-1,        6,8,-1,        8,8,-2,
20765   DATA 9,8,-1,        11,8,-2,       16,8,-1,       18,8,-2,
20770   DATA 19,8,-1,       21,8,-2,       27,8,-1
20775  Psil: DATA 1,12,-2,   1,12,-1,       2,11,-1,       2,8,-1,
20780   DATA 3,7,-1,      5,7,-1,        6,8,-1,        6,11,-1,
20785   DATA 7,12,-1,     3,12,-2,       3,12,-1,       5,12,-1,
20790   DATA 4,12,-2,     4,12,-1,       4,4,-1,        3,4,-2,
20795   DATA 3,4,-1,      5,4,-1
20800  Phi: DATA 1,6,-2,     1,6,-1,        3,5,-1,        5,5,-1,
20805   DATA 7,6,-1,      7,8,-1,        5,9,-1,        3,9,-1,
20810   DATA 1,8,-1,      1,6,-1,        4,11.5,-2,     4,11.5,-1,
20815   DATA 4,2.5,-1
20820  Tau: DATA 1,8,-2,     1,8,-1,        2,9,-1,        6,9,-1,
20825   DATA 4,9,-2,      4,9,-1,        4,4,-1,        6,4,-1
20830  Theta: DATA 1,7.5,-2, 1,7.5,-1,      5,7.5,-1,      5,9,-1,
20835   DATA 4,11,-1,     2,11,-1,       1,9,-1,        1,6,-1,
20840   DATA 2,4,-1,      4,4,-1,        5,6,-1,        5,7.5,-1
20845  Lamda: DATA 1,4,-2,  3.7,8,-1,      1,11,-2,       2,11,-1,
20850   DATA 6,4,-1,      7,4,-1
20855  Omegal: DATA 1,5,-2,   1,5,-1,        1,4,-1,        3,4,-1,
20860   DATA 3,6,-1,      1,8,-1,        1,10,-1,       3,12,-1,
20865   DATA 5,12,-1,     7,10,-1,       7,8,-1,        5,6,-1,
20870   DATA 5,4,-1,      7,4,-1,        7,5,-1
20875  Eta: DATA 1,10,-2,    1,10,-1,       2,9,-1,        2,4,-2,
20880   DATA 2,4,-1,      2,9,-1,        3,10,-1,       5,10,-1,
20885   DATA 6,8,-1,      6,2,-1
20890  Pai: DATA 1,9,-2,     1,9,-1,        2,10,-1,       7,10,-1,
20895   DATA 3,10,-2,     3,10,-1,       3,4,-1,        6,10,-2,
20900   DATA 6,10,-1,     6,4,-1
20905  Intgrl: DATA 2,0,-2,  2,0,-1,        3,0,-1,        4,1,-1,
20910   DATA 5,14,-1,     6,15,-1,       7,15,-1
20915  Sbs: DATA 2,3,-2,    2,8,-1,        2,5,-2,        3,6,-1
20920   DATA 4,6,-1,    5,5,-1,        5,4,-1,        4,3,-1
20925   DATA 3,3,-1,    2,4,-1
20930  S0: DATA 2,4,-2,     2,6,-1,        3,7,-1,        4,7,-1
20935   DATA 5,6,-1,     5,4,-1,        4,3,-1,        3,3,-1
20940   DATA 2,4,-1
20945  S1: DATA 2,6,-2,    3,7,-1,        3,3,-1,        2,3,-2
20950   DATA 4,3,-1
20955  Circle: DATA 2,7,-2,  2,8,-1,      2.7,9.2,-1,    4,10,-1
20960   DATA 5,10,-1,      6.3,9.2,-1,  7,8,-1,        7,7,-1
20965   DATA 6.3,5.8,-1,   5,5,-1,      4,5,-1,        2.7,5.8,-1
20970   DATA 2,7,-1
20975  Cross: DATA 1,4,-2,   8,11,-1,     1,11,-2,       8,4,-1
20980  Tangl: DATA 1,5,-2,   4.5,11,-1,   8,5,-1,        1,5,-1
20985  Square: DATA 2,5,-2,  2,10,-1,     7,10,-1,       7,5,-1
20990   DATA 2,5,-1
20995  Dot: DATA 2,7,-2,  2,8,-1,      2.7,9.2,-1,    4,10,-1
21000   DATA 5,10,-1,      6.3,9.2,-1,  7,8,-1,        7,7,-1
21005   DATA 6.3,5.8,-1,   5,5,-1,      4,5,-1,        2.7,5.8,-1
21010   DATA 2,7,-1,       3,9,-1,      3,6,-1,        4,10,-1
21015   DATA 4,5,-1,       5,10,-1,     5,5,-1,        6,9,-1
21020   DATA 6,6,-1
21025   ! ------------------Replace-----------------------------
21030   Old_chars$=""     ! In case anything is left in COM from the last run...
21035   New_udc(CHR$(168),Sigma(*))
21040   New_udc(CHR$(169),Infinity(*))
21045   New_udc(CHR$(170),Arrowb(*))
21050   New_udc(CHR$(171),Box(*))
21055   New_udc(CHR$(172),Alpha(*))
21060   New_udc(CHR$(174),Deg(*))
21065   New_udc(CHR$(173),Degc(*))
21070   New_udc(CHR$(224),Permill(*))
21075   New_udc(CHR$(228),Delta(*))
21080   New_udc(CHR$(225),Solid(*))
21085   New_udc(CHR$(226),Dotted(*))
21090   New_udc(CHR$(227),Chain(*))
21095   New_udc(CHR$(229),Psil(*))
21100   New_udc(CHR$(230),Phi(*))
21105   New_udc(CHR$(150),Rho(*))
21110   New_udc(CHR$(231),Tau(*))
21115   New_udc(CHR$(232),Theta(*))
21120   New_udc(CHR$(151),Lamda(*))
21125   New_udc(CHR$(233),Omegal(*))
21130   New_udc(CHR$(234),Pai(*))
21135   New_udc(CHR$(235),Eta(*))
21140   New_udc(CHR$(236),Intgrl(*))
21145   New_udc(CHR$(237),Arrow3(*))
21150   New_udc(CHR$(238),Sbs(*))
21155   New_udc(CHR$(240),S0(*))
21160   New_udc(CHR$(241),S1(*))
21165   New_udc(CHR$(239),Arrow8(*))
21170   New_udc(CHR$(180),Circle(*))
21175   New_udc(CHR$(181),Dot(*))
21180   New_udc(CHR$(182),Cross(*))
21185   New_udc(CHR$(183),Tangl(*))
21190   New_udc(CHR$(184),Square(*))
21195   GOSUB Menu
21200   SUBEXIT
21205  Menu:! ------------------Menu--------------------------------
21210   OUTPUT 1 USING "/"
21215   OUTPUT 1;"******************* Menu of symbol ********************"
21220   OUTPUT 1;"GRE_S: Alpha(172),Phi(230),Tau(231),Theta(232),Eta(235)"
21225   OUTPUT 1;"       Rho(150),Lamda(151)"
21230   OUTPUT 1;"GRE_L: Psil(229),Omegal(233)"
21235   OUTPUT 1;"MATHE: Sigmal(168),Pai(234),Deltal(228),Integral(236),Infinity(169)"
21240   OUTPUT 1;"FIGUR: Arrowb(170),Arrow3(237),Arrow8(239),Box(171)"
21245   OUTPUT 1;"UNIT : Deg(174),Degc(173),Permil(224)"
21250   OUTPUT 1;"LINE : Solid(225),Dotted(226),Chain(227)"
21255   OUTPUT 1;"SUFAP: Sbs(238)"
21260   OUTPUT 1;"SUFNM: S0(240),S1(241)"
21265   OUTPUT 1;"SYMBL: Circle(180),Dot(181),Cross(182),Tangl(183),Square(184)"
21270   OUTPUT 1 USING "/"
21275   RETURN
21280 SUBEND
21285 ! ####################################################
21290 SUB Seqrpfile
21295   ! Preparation for reading of plot file made by Pltfile.
21300   ! ####################################################
21305   COM /Seqrpf/Basefname$[8],Fnmin,Fnmax,Op$[8],Pr$[8],Py$[8],Pl$[8],Fnum,Endflg   ! -->
21310   PRINT
21315   PRINT "Sequencial reading of plot file."
21320   PRINT "********************************"
21325   DISP "Input basic file name."
21330   WAIT 2
21335   Filename(Basefname$)              ! Input basic file name.
21340   PRINT "Basic file name:",Basefname$
21345   INPUT "File#: min,max (00-99) ?",Fnmin$,Fnmax$
21350   PRINT "File#: min,max ",Fnmin$,Fnmax$
21355   Fnmin=VAL(Fnmin$)
21360   Fnmax=VAL(Fnmax$)
21365   INPUT "Print para. ? (Y/N)",Op$
21370   PRINT "Print para.",Op$
21375   INPUT "Print Y data (Y/N) ?",Py$
21380   PRINT "Print Y data.",Py$
21385   INPUT "Printer (C/P) ?",Pr$
21390   PRINT "Printer (C/P)",Pr$
21395   INPUT "Plot (Y/N) ?",Pl$
21400   PRINT "Plot ",Pl$
21405   IF Pl$="Y" THEN CALL Seqplot
21410   INPUT "Input OK (Y/N) ?",Ok$
21415   IF Ok$="N" THEN GOTO 21335
21420   Fnum=0                            ! Reset. Counted in Pltfile.
21425   Endflg=0                          ! Reset. Senced in Pltfile.
21430 SUBEND
21435 ! ####################################################
21440 SUB Tsavef
21445 ! NOF averaging of theoretical TS functions calculated & filed by
21450 !   Seqbackpat and read by Seqrpfile & Pltfile.
21455 ! Can superpose measured data by Plotpoint(soft) or Any_posi_led(liquid).
21460 ! ####################################################
21465   COM /Pltd/Lnum,Dnum,Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,Auto,Autled$,Linlbl$(*),Manuled   ! <-- Pltfile
21470   COM /Seqrpf/Basefname$[8],Fnmin,Fnmax,Op$[8],Pr$[8],Py$[8],Pl$[8],Fnum,Endflg   ! <-- Seqrpfile
21475   ! ------------------Input-------------------------------
21480   DIM Ori(1:2,1:10)         ! Mean_S.D., #
21485   PRINT "NOF averaging of theoretical TS functions."
21490   PRINT "******************************************"
21495   INPUT "Superpose measured data when soft (Y/N) ?",Sm$
21500   IF Sm$="Y" THEN
21505     INPUT "Block(Speicies) number ?",Bn
21510     ALLOCATE Symblb$(1:Bn)[20]
21515     FOR I=1 TO Bn
21520       DISP "Symbol label #";I,
21525       INPUT Symblb$(I)
21530     NEXT I
21535   END IF
21540  Oori: PRINT
21545   PRINT "Orientation","#","Mean","S.D."
21550   Io=1
21555  I: INPUT "Orientation dist.: Mean, S.D. (1_dir.:*,0  End:*,-1) ?",Meant,St
21560   IF St=-1 THEN
21565     Iomax=Io-1
21570     GOTO 21610
21575   ELSE
21580     PRINT " "," ",Io,Meant,St
21585     Ori(1,Io)=Meant
21590     Ori(2,Io)=St
21595     Io=Io+1
21600     GOTO I
21605   END IF
21610   PRINT
21615   INPUT "Bladder inc. angle (Ex.6) ?",Bang
21620   PRINT "Bladder inc. angle ",Bang
21625   ! ------------------File in & Ave.----------------------
21630   Seqrpfile
21635   Ihmax=Fnmax-Fnmin+1
21640   Ih=1
21645   ALLOCATE Tsa(1:2,1:Iomax,1:Ihmax)      ! Kind,Orient.,H
21650  H: Pltfile(1)
21655   IF Endflg=1 THEN GOTO Res        ! End of input & ave.
21660   Icent=INT(Dnum/2)+1              ! Center angle index.
21665   Ib=INT(Bang/Xinc+.1)
21670   ALLOCATE Ts(1:Dnum)              ! For FNTsnof.
21675   FOR K=1 TO 2                     ! Soft,Liquid
21680     FOR It=1 TO Dnum
21685       Ts(It)=Y(K-1,It-1)           ! Option base of Pltfile = 0
21690     NEXT It
21695     FOR Io=1 TO Iomax
21700       IF Ori(2,Io)=0 THEN          ! 1 direction
21705         Iang=INT(Ori(1,Io)/Xinc+.1)
21710         Tsa(K,Io,Ih)=Ts(Icent+Iang+(K=1)*Ib)
21715       ELSE
21720         Tsa(K,Io,Ih)=FNTsnof(Ts(*),Xinc,Ori(1,Io)+Bang*(K=1),Ori(2,Io))
21725         ! Bladder inclination is realized by shift of Meant.
21730       END IF
21735     NEXT Io
21740   NEXT K
21745   DEALLOCATE Ts(*)
21750   Ih=Ih+1
21755   GOTO H
21760  Res:   ! ------------------Print-------------------------------
21765   INPUT "H0: Min, Max, Step ?",H0min,H0max,H0step
21770   INPUT "b/a, ab/a ?",Ba,Aba
21775   Kl=PI*SQR(1-Ba^2)              ! PI/Xi
21780   INPUT "Rho1/0, C1/0 ?",R10,C10
21785   FOR K=1 TO 2
21790     PRINT
21795     IF K=1 THEN PRINT "Soft"
21800     IF K=2 THEN PRINT "Liquid"
21805     PRINT " File#",
21810     FOR Io=1 TO Iomax
21815       PRINT "(";Ori(1,Io);",";Ori(2,Io);")",
21820     NEXT Io
21825     PRINT
21830     Ih=1
21835     FOR Fnum=Fnmin TO Fnmax
21840       PRINT Fnum,
21845       FOR Io=1 TO Iomax
21850         Tsa(K,Io,Ih)=Tsa(K,Io,Ih)-40+(K=1)*20*LGT(Aba)   ! -40 to cm.
21855         PRINT PROUND(Tsa(K,Io,Ih),-2),
21860       NEXT Io
21865       PRINT
21870       Ih=Ih+1
21875     NEXT Fnum
21880   NEXT K
21885   ! ------------------Interp. & Plot----------------------
21890   DIM Dom(1:200),Al$[120]
21895   ALLOCATE Xs(1:Ihmax),Ys(1:Ihmax)    ! For Spline
21900   INPUT "H0 step to be interpolated ?",H0stepi
21905   INPUT "Linear or Semi-log (LIN/SML) ?",G$
21910   Ihi=1
21915   FOR H0=H0min TO H0max STEP H0stepi
21920     Dom(Ihi)=H0
21925     Ihi=Ihi+1
21930     H0maxi=H0
21935   NEXT H0
21940   Narg=Ihi-1
21945   REDIM Dom(1:Narg)
21950   ALLOCATE Fun(1:Narg),Der(1:Narg),Yp(1:Iomax,1:Narg),Ll$(1:Iomax)[20]
21955   H0=H0min
21960   FOR Ih=1 TO Ihmax
21965     Xs(Ih)=H0
21970     H0=H0+H0step
21975   NEXT Ih
21980   Al$="Tsavef "&Basefname$&VAL$(Fnmin)&"-"&VAL$(Fnmax)&" b/a="&VAL$(Ba)
21985   FOR K=1 TO 2
21990     IF K=1 THEN
21995       Al$=Al$&" Soft"&" ab/a="&VAL$(Aba)&" Bang="&VAL$(Bang)
22000       Fhl=Kl*Aba                 ! PI/Xi*ab/a
22005     ELSE
22010       Al$=Al$&" Liquid"&" R10="&VAL$(R10)&" C10="&VAL$(C10)
22015       Fhl=Kl
22020     END IF
22025     FOR Io=1 TO Iomax
22030       FOR Ih=1 TO Ihmax
22035         Ys(Ih)=Tsa(K,Io,Ih)
22040       NEXT Ih
22045       Spline(Ihmax,Narg,Xs(*),Ys(*),Dom(*),Fun(*),Der(*),Int,1.E-4)
22050       FOR Ihi=1 TO Narg
22055         Yp(Io,Ihi)=Fun(Ihi)
22060       NEXT Ihi
22065       Ll$(Io)="("&VAL$(Ori(1,Io))&","&VAL$(Ori(2,Io))&")"
22070       IF (K=1 AND Ori(1,Io)=-Bang AND Ori(2,Io)=0) OR (K=2 AND Ori(1,Io)=0 AND Ori(2,Io)=0) THEN Ll$(Io)=Ll$(Io)&" MAX."
22075     NEXT Io
22080     Symbol(1)
22085     Hstep=H0stepi/Fhl
22090     IF G$="SML" THEN Hstep=-Hstep
22095     Plot(Yp(*),H0min/Fhl,H0maxi/Fhl,Hstep,G$,"L/"&CHR$(151)," A dB"," ",0,Al$,Ll$(*),(K=2))
22100     IF Sm$="Y" AND K=1 THEN
22105       DISP "Change file to measured data"
22110       WAIT 2
22115       Plotpoint(G$," "," ","SPP"," ",Symblb$(*),0)
22120     END IF
22125     IF K=1 THEN
22130       DISP "Next for liquid."
22135       WAIT 2
22140     END IF
22145   NEXT K
22150   INPUT "Any other orientation ? (Y/N)",Aop$
22155   IF Aop$="Y" THEN GOTO Oori
22160 SUBEND
22165 ! ####################################################
22170 DEF FNGauss(X,M,S)
22175 ! Normal probability density.
22180 ! ####################################################
22185   RETURN EXP(-(X-M)*(X-M)/2/S/S)/S/SQR(2*PI)
22190 FNEND
22195 ! ####################################################
22200 DEF FNSimp(Y(*),H)
22205 ! Numerical integration by Simpson's rule.
22210 ! Data are passed by array.
22215 ! ####################################################
22220   Imin=BASE(Y,1)
22225   Size=SIZE(Y,1)
22230   IF Size<3 THEN
22235     BEEP
22240     PRINT "Parameter error in FNSimp. Size<3. Return 0."
22245     RETURN 0
22250   END IF
22255   IF Size MOD 2=0 THEN
22260     BEEP
22265     PRINT "Warning in FNSimp. Size is even. Last data is not used."
22270     Size=Size-1
22275   END IF
22280   Imax=Imin+Size-1
22285   Int=0
22290   FOR I=Imin TO Imax-2 STEP 2
22295     Int=Int+Y(I)+4*Y(I+1)+Y(I+2)
22300   NEXT I
22305   RETURN Int*H/3
22310 FNEND
22315 ! ####################################################
22320 SUB Spline(N,Narg,X(*),Y(*),Domain(*),Func(*),Deriv(*),Int,Eps)
22325 ! Interpolation  by Spline method.  See math. lib. of HP.
22330 ! ####################################################
22335   OPTION BASE 1
22340   ALLOCATE S(N),G(N-1),Work(N-1)
22345   FOR I=2 TO N-1
22350     Xi=X(I)
22355     Xim1=X(I-1)
22360     Xip1=X(I+1)
22365     Yi=Y(I)
22370     Yim1=Y(I-1)
22375     Yip1=Y(I+1)
22380     Xd=Xi-Xim1
22385     H=Xip1-Xim1
22390     Work(I)=.5*Xd/H
22395     T=((Yip1-Yi)/(Xip1-Xi)-(Yi-Yim1)/Xd)/H
22400     S(I)=2*T
22405     G(I)=3*T
22410   NEXT I
22415   S(1)=0
22420   S(N)=0
22425   W=8-4*SQR(3)
22430   U=0
22435   FOR I=2 TO N-1
22440     T=W*(-S(I)-Work(I)*S(I-1)-(.5-Work(I))*S(I+1)+G(I))
22445     H=ABS(T)
22450     IF H>U THEN U=H
22455     S(I)=S(I)+T
22460   NEXT I
22465   IF U>=Eps THEN 22430
22470   FOR I=1 TO N-1
22475     G(I)=(S(I+1)-S(I))/(X(I+1)-X(I))
22480   NEXT I
22485   IF Narg=0 THEN 22595
22490   FOR J=1 TO Narg
22495  Corrector: I=1
22500     T=Domain(J)
22505     IF T>=X(1) THEN 22525
22510     BEEP
22515     PRINT "Error in Spline. Argument out of bounds."
22520     SUBEXIT
22525     I=I+1
22530     IF I>N THEN 22510
22535     IF T>X(I) THEN 22525
22540     I=I-1
22545     H=Domain(J)-X(I)
22550     T=Domain(J)-X(I+1)
22555     Xd=H*T
22560     Sp=S(I)+H*G(I)
22565     Z=1/6
22570     U=Z*(S(I)+S(I+1)+Sp)
22575     W=(Y(I+1)-Y(I))/(X(I+1)-X(I))
22580     Func(J)=W*H+Y(I)+Xd*U
22585     Deriv(J)=W+(H+T)*U+Z*Xd*G(I)
22590   NEXT J
22595   Int=0
22600   FOR I=1 TO N-1
22605     H=X(I+1)-X(I)
22610     Int=Int+.5*H*(Y(I)+Y(I+1))-1-24*H^3*(S(I)+S(I+1))
22615   NEXT I
22620 SUBEND
22625 ! ####################################################
22630 DEF FNTsnof(Tsfun(*),Tinc,Meant,Sigmat)
22635 ! Nakken-Olsen-Foote averaging of TS.
22640 ! Tsfun is TS in dB as function of pitch angle symmetric with 0 deg.
22645 ! Results in dB.
22650 ! Orientation distribution is Gaussian with Meant and Sigmat truncated
22655 !  at Sigmat*3
22660 ! If Tsfun is theoretical for soft spheroid , Meant must be
22665 !  True_meant+ABS(Bladder_inc_ang).
22670 ! ####################################################
22675   Jb=BASE(Tsfun,1)
22680   Js=SIZE(Tsfun,1)
22685   Jmax=Jb+Js-1
22690   Tr=INT(Js/2+.1)*Tinc       ! Tsfun is given from -Tr TO Tr.
22695   Tg=Sigmat*3                ! Gaussian pdf is truncated at Tg.
22700   T1=Meant-Tg
22705   T2=Meant+Tg
22710   Imax=INT(2*Tg/Tinc+.1)+1
22715   ALLOCATE Y(1:Imax)
22720   IF T1<-Tr THEN             ! Average of 1st 10deg
22725     Sum=0
22730     I=0
22735     T=0
22740     WHILE T<=10
22745       Sum=Sum+10^(Tsfun(Jb+I)/10)
22750       T=T+Tinc
22755       I=I+1
22760     END WHILE
22765     Tslow=Sum/(I-1)
22770   END IF
22775   IF T2>Tr THEN             ! Average of last 10deg
22780     Sum=0
22785     I=0
22790     T=0
22795     WHILE T<=10
22800       Sum=Sum+10^(Tsfun(Jmax-I)/10)
22805       T=T+Tinc
22810       I=I+1
22815     END WHILE
22820     Thigh=Sum/(I-1)
22825   END IF
22830   I=0
22835   FOR T=T1 TO T2 STEP Tinc
22840     I=I+1
22845     SELECT T
22850     CASE <-Tr
22855       Ts=Tslow
22860     CASE >Tr
22865       Ts=Tshigh
22870     CASE ELSE
22875       J=Jb+INT((T+Tr)/Tinc+.1)
22880       Ts=10^(Tsfun(J)/10)
22885     END SELECT
22890     Y(I)=Ts*FNGauss(T,Meant,Sigmat)
22895   NEXT T
22900   RETURN 10*LGT(FNSimp(Y(*),Tinc)/.997)
22905 FNEND
22910 ! ####################################################
22915 SUB Drawmove(G$,Xx,Yy,Dm,OPTIONAL Linetype,Rep)
22920 ! Draw(Dm=0) or Move(1) to Xx,Yy.
22925 ! ####################################################
22930   X=Xx
22935   Y=Yy
22940   IF G$="SML" OR G$="LOG" THEN X=LGT(X)
22945   IF G$="LOG" THEN Y=LGT(Y)
22950   SELECT Dm
22955   CASE 0
22960     LINE TYPE Linetype,Rep
22965     DRAW X,Y
22970     LINE TYPE 1
22975   CASE 1
22980     MOVE X,Y
22985   END SELECT
22990 SUBEND
22995 ! ####################################################
23000 SUB Plotpoint(G$,Xlb$,Ylb$,Ttl$,Autled$,Symblb$(*),Manuled)
23005 !  Plotting of points filed by blocking mode of Filein.
23010 !  Format in Filein
23015 !    Blocking. Block# is kind#.
23020 !    There must be only one series of read data, i.e. (Item,Type)=(QES,R).
23025 !    Odd index is X data, even Y.
23030 !    Para(1)=Stopgap>Xmax.
23035 !  Can superpose on graph by Plot through Ttl$="SPP".
23040 ! ####################################################
23045   OPTION BASE 1
23050   DEG
23055   COM /Pltr/Plot$[8],Sp$[8],Lastpen   ! Plot$<--Graphinit
23060   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl  ! <-- Axlas
23065   COM /Lgnd/Le(1:5,1:20,1:4),Legen$(1:5,1:20)[20],Lnn(1:5),Fn,Fnmax  ! <-- Any_posi_led
23070   COM /Gdu/Xgmax,Ygmax,Xupg,Yupg    ! Upg<--, Max-->
23075   DIM Filename$[10],Com$[100],Para$(20)[20],Para(20),Item$(500,2)[10],Rd(10,500),Cd$(1,1)[10]
23080   INTEGER Id(1,1),Nk,Np,Nr,Ni,Nc
23085   Cslsl=Csl*.7                      ! Used for symbol label
23090   Cslsb=Csl*.5                      ! Make symbol small.
23095   Symblflg=0
23100   ! -----------------------Input--------------------------
23105   Fileio("R",Filename$,Com$,Nk,Np,Nr,Ni,Nc,Para$(*),Para(*),Item$(*),Rd(*),Id(*),Cd$(*))
23110   IF Nk>5 THEN Symblflg=1
23115   IF Nk>8 THEN
23120     BEEP
23125     PRINT "Too many symbol."
23130   END IF
23135   Imax=INT(Nr/2+.1)
23140   ALLOCATE Xy(Nk,Imax,2)
23145   FOR K=1 TO Nk
23150     FOR Ir=1 TO Nr-1 STEP 2
23155       I=INT(Ir/2)+1
23160       Xy(K,I,1)=Rd(K,Ir)           ! X data
23165       Xy(K,I,2)=Rd(K,Ir+1)         ! Y data
23170     NEXT Ir
23175   NEXT K
23180   Stopgap=0  !*********
23185   ! ------------------Processing--------------------------
23190   INPUT "Print(P),Correct(C),Modify(M),No(N) ?",Pcmn$
23195  Fmt: IMAGE "K ",DD,3X,"I ",DDD,3X,"X ",SD.DDDE,3X,"Y ",SD.DDDE
23200   SELECT Pcmn$
23205   CASE "P"
23210     PRINT
23215     PRINT "X,Y data"
23220     PRINT "********"
23225     PRINT
23230     FOR K=1 TO Nk
23235       FOR I=1 TO Imax
23240         PRINT USING Fmt;K,I,Xy(K,I,1),Xy(K,I,2)
23245       NEXT I
23250       PRINT
23255     NEXT K
23260   CASE "C"
23265     PRINT
23270     INPUT "Which data (K,I) ?",K,I
23275     PRINT USING Fmt;K,I,Xy(K,I,1),Xy(K,I,2)
23280     INPUT "New data (X,Y) ?",Xy(K,I,1),Xy(K,I,2)
23285     PRINT USING Fmt;K,I,Xy(K,I,1),Xy(K,I,2)
23290     INPUT "More correction (Y/N) ?",Mc$
23295     IF Mc$="Y" THEN GOTO 23265
23300   CASE "M"
23305     INPUT "Xp,Xm,Yp,Ym of (Xdata+Xp)*Xm & (Ydata+Yp)*Ym ?",Xp,Xm,Yp,Ym
23310     FOR K=1 TO Nk
23315       FOR I=1 TO Imax
23320         Xy(K,I,1)=(Xy(K,I,1)+Xp)*Xm
23325         Xy(K,I,2)=(Xy(K,I,2)+Yp)*Ym
23330       NEXT I
23335     NEXT K
23340   CASE "N"
23345     GOTO 23360
23350   END SELECT
23355   GOTO 23190
23360   READ Xmin,Xmax,Ymin,Ymax
23365   DATA 1E30,-1E30,1E30,-1E30
23370   FOR K=1 TO Nk
23375     FOR I=1 TO Imax
23380       IF Xy(K,I,1)=Stopgap THEN Nextk
23385       Xmin=MIN(Xy(K,I,1),Xmin)
23390       Xmax=MAX(Xy(K,I,1),Xmax)
23395       Ymin=MIN(Xy(K,I,2),Ymin)
23400       Ymax=MAX(Xy(K,I,2),Ymax)
23405     NEXT I
23410  Nextk: NEXT K
23415   ! ------------------Mode & Scaling----------------------
23420   INPUT "Only points(P)/ with polygonal line(L)/ with reg. line(R) ?",Pl$
23425   IF Pl$="L" THEN
23430     ALLOCATE Lt(Nk,2)
23435  Lt: DATA 1,0,4,1,6,5,5,3,4,2,8,8,7,8,4,4,6,5,5,2,8,5,7,5,4,3,6,3,5,5                                         ! For line type.
23440     INPUT "Line type: Bold(B)/Dashed(D)/Classify(C) ?",Bdc$
23445     IF Bdc$="C" THEN
23450       RESTORE Lt
23455       FOR K=1 TO Nk
23460         MAT SORT Xy(K,*,1)
23465       NEXT K
23470     END IF
23475     FOR K=1 TO Nk
23480       SELECT Bdc$
23485       CASE "B"
23490         Lt(K,1)=1
23495         Lt(K,2)=0
23500       CASE "D"
23505         Lt(K,1)=5
23510         Lt(K,2)=3
23515       CASE "C"
23520         READ Lt(K,1),Lt(K,2)
23525       CASE ELSE
23530         GOTO 23440
23535       END SELECT
23540     NEXT K
23545   END IF
23550   IF Ttl$="SPP" THEN 23560              ! Superpose on another scaling
23555   Axlas(G$,Xmin,Xmax,Ymin,Ymax,1,Xlb$,Ylb$,Ttl$,0)
23560   ! ------------------Plot points-------------------------
23565   Symbol(1)
23570   RESTORE Pp
23575  Pp: DATA 180,181,182,183,184
23580   PEN 1
23585   FOR K=1 TO Nk
23590     GOSUB Readsymb
23595     FOR I=1 TO Imax
23600       X=Xy(K,I,1)
23605       Y=Xy(K,I,2)
23610       !       PRINT "X,Y",X,Y  !******
23615       IF X=Stopgap THEN GOTO 23650
23620       IF Pl$="L" AND I<>1 THEN CALL Drawmove(G$,X,Y,0,Lt(K,1),Lt(K,2))
23625       LORG 5
23630       Drawmove(G$,X,Y,1)
23635       Label(CHR$(Symb))
23640       IF Pl$="L" THEN CALL Drawmove(G$,X,Y,1)
23645     NEXT I
23650   NEXT K
23655   LDIR 0
23660      ! ------------------Symbol label------------------------
23665   RESTORE Pp
23670   Ycs=Cslsl*Yupg
23675   Xcs=Cslsl*9/15*Xupg
23680  Otherposi: INPUT "Position of symbol label(UL,UR,BL,BR,NO) ?",Posi$
23685   IF Plot$="705" AND Posi$<>"NO" THEN INPUT "Plot symbol name(Y/N) ?",Pln$
23690   IF Posi$="BR" OR Posi$="UR" THEN
23695     Maxlen=0
23700     FOR K=1 TO Nk
23705       IF LEN(Symblb$(K))>Maxlen THEN Maxlen=LEN(Symblb$(K))
23710     NEXT K
23715     Len=Maxlen+3
23720   END IF
23725   SELECT Posi$
23730   CASE "UL"
23735     Xsp=Xsmin+Xsr/25
23740     Ysp=Ysmax-Ysr/15
23745   CASE "UR"
23750     Xsp=Xsmax-Xcs*Len
23755     Ysp=Ysmax-Ysr/15
23760   CASE "BL"
23765     Xsp=Xsmin+Xsr/25
23770     Ysp=Ysmin+Ycs*Nk
23775   CASE "BR"
23780     Xsp=Xsmax-Xcs*Len
23785     Ysp=Ysmin+Ycs*Nk
23790   CASE "NO"
23795     GOTO 23900
23800   CASE ELSE
23805     GOTO 23680
23810   END SELECT
23815   Yl=Ysp+Ycs
23820   FOR K=1 TO Nk
23825     GOSUB Readsymb
23830     Yl=Yl-Ycs
23835     MOVE Xsp,Yl
23840     PEN 1
23845     LORG 5
23850     Label(CHR$(Symb))
23855     MOVE Xsp+1*Xcs,Yl
23860     PEN 2
23865     LORG 2
23870     CSIZE Cslsl
23875     LDIR 0
23880     Label(Symblb$(K))
23885   NEXT K
23890   LDIR 0
23895   ! ------------------Legend------------------------------
23900   PEN 2
23905   IF Autled$=" " THEN GOTO 23925
23910   Wl$="Y"
23915   IF Plot$="705" THEN INPUT "Need auto legend (Y/N) ?",Wl$
23920   IF Wl$="Y" THEN CALL Autolegend(Autled$)
23925   IF Manuled=1 THEN
23930     Fn=0
23935     Legend
23940   END IF
23945   IF Plot$="701" THEN DUMP GRAPHICS #701
23950   ! ------------------For next work-----------------------
23955   INPUT "Other plot or exit (P/E) ?",Pe$
23960   IF Pe$="P" THEN
23965     RESTORE Pp
23970     GCLEAR
23975     GOTO 23420
23980   ELSE
23985     SUBEXIT
23990   END IF
23995   SUBEXIT
24000  Readsymb:! ------------------------------------------------------
24005   IF K=6 THEN
24010     RESTORE Pp
24015     READ Symb                   ! Pass circle
24020     READ Symb                   !      dot
24025   END IF
24030   IF K>=6 THEN LDIR 45          ! Make other symbols
24035   READ Symb
24040   PRINT "Symb,1025",Symb
24045   CSIZE Cslsb
24050   IF K=2 THEN CSIZE Cslsb*.5
24055   RETURN
24060 SUBEND
24065 ! ####################################################
24070 SUB Filein
24075 !  Input blocked or unblocked data and write to and read from file.
24080 !  Unblocked data are only one type, therefore use Para properly.
24085 !  If only one block in blocked mode, block 2 must be dummy.
24090 ! ####################################################
24095   OPTION BASE 1
24100   INTEGER Nb,Np,Nr,Ni,Nc
24105   INPUT "Printer(C/P) ?",Prt$
24110   IF Prt$="P" THEN PRINTER IS 701
24115   DIM Comment$[100]
24120   LINPUT "Comment ?",Comment$
24125   INPUT "Need Parameter, Blocking(Y/N) ?",Para$,Block$
24130   PRINT
24135   PRINT "**************************************************"
24140   PRINT "Comment"
24145   PRINT TAB(2),Comment$
24150   ! ------------------Parameter---------------------------
24155   IF Para$="N" THEN
24160     ALLOCATE Paraname$(1)[20],Paraval(1)
24165     Np=0
24170   ELSE
24175     ALLOCATE Paraname$(20)[20],Paraval(20)
24180     Ip=1
24185  Para: DISP "Parameter ";Ip;": Name, Value(<TSAL,0> to end)",
24190     INPUT Paraname$(Ip),Paraval(Ip)
24195     IF Paraname$(Ip)="TSAL" THEN
24200       Np=Ip-1
24205       REDIM Paraname$(Np),Paraval(Np)
24210     ELSE
24215       Ip=Ip+1
24220       GOTO Para
24225     END IF
24230     PRINT
24235     PRINT "Parameters"
24240     FOR Ip=1 TO Np
24245       PRINT TAB(2),Paraname$(Ip),Paraval(Ip)
24250     NEXT Ip
24255   END IF
24260 ! ------------------Data name & type--------------------
24265   IF Block$="N" THEN
24270     INPUT "Approx. max. data num. ?",Dnmax
24275     Nb=1
24280     ALLOCATE Item$(1,2)[10]
24285     INPUT "Data Name, Type(R/I/C) ?",Item$(1,1),Item$(1,2)
24290     SELECT Item$(1,2)
24295     CASE "R"
24300       ALLOCATE REAL Rdata(1,Dnmax)
24305       ALLOCATE INTEGER Idata(1,1)     ! Dummy
24310       ALLOCATE Cdata$(1,1)[10]        ! Dummy
24315     CASE "I"
24320       ALLOCATE INTEGER Idata(1,Dnmax)
24325       ALLOCATE REAL Rdata(1,1)
24330       ALLOCATE Cdata$(1,1)[10]
24335     CASE "C"
24340       ALLOCATE Cdata$(1,Dnmax)[10]
24345       ALLOCATE REAL Rdata(1,1)
24350       ALLOCATE INTEGER Idata(1,1)
24355     END SELECT
24360     PRINT
24365     PRINT "Data name & Type",Item$(1,1),Item$(1,2)
24370   ELSE
24375     INPUT "Approx. max. num.:Block,Item ?",Bnmax,Inmax
24380     ALLOCATE Item$(Inmax,2)[10]          ! Item name, Type
24385     ALLOCATE INTEGER Itable(Inmax)           ! Item # to # each type
24390     ALLOCATE Seq(Inmax,2)
24395     DIM Dispseq$[30]
24400     MAT Seq=(0)
24405     Id=1
24410     Ir=0
24415     Ii=0
24420     Ic=0
24425  Nextitem: DISP "Item ";Id;": Name, Type(R/I/C) (<TSAL,R> to end, <QES,Type> to seq.)",
24430     INPUT Item$(Id,1),Item$(Id,2)
24435     SELECT Item$(Id,1)
24440     CASE "TSAL"
24445       Nd=Id-1
24450       Nr=Ir
24455       Ni=Ii
24460       Nc=Ic
24465       REDIM Item$(Nd,2),Itable(Nd)
24470     CASE "QES"
24475       Type$=Item$(Id,2)
24480       INPUT "Name, Start#, End#, Stopgap(Ex.-32768) ?",Name$,N1,N2,Sg
24485       FOR Num=N1 TO N2
24490         Item$(Id,1)=Name$&VAL$(Num)
24495         Item$(Id,2)=Type$
24500         Seq(Id,1)=1                    ! Stopgap flag for each item.
24505         Seq(Id,2)=Sg
24510         GOSUB Inc
24515       NEXT Num
24520       GOTO Nextitem
24525     CASE ELSE
24530       GOSUB Inc
24535       GOTO Nextitem
24540     END SELECT
24545     PRINT
24550     PRINT "Data name & Type"
24555     PRINT "No.","Name","Type","No. in type"
24560     FOR Id=1 TO Nd
24565       PRINT Id,Item$(Id,1),Item$(Id,2),Itable(Id)
24570     NEXT Id
24575     IF Nr<>0 THEN
24580       ALLOCATE REAL Rdata(Bnmax,Nr)
24585     ELSE
24590       ALLOCATE REAL Rdata(1,1)
24595     END IF
24600     IF Ni<>0 THEN
24605       ALLOCATE INTEGER Idata(Bnmax,Ni)
24610     ELSE
24615       ALLOCATE INTEGER Idata(1,1)
24620     END IF
24625     IF Nc<>0 THEN
24630       ALLOCATE Cdata$(Bnmax,Nc)[10]
24635     ELSE
24640       ALLOCATE Cdata$(1,1)[10]
24645     END IF
24650   END IF
24655 ! ------------------Input blocked data------------------
24660   Ib=1
24665   PRINT
24670   PRINT "Data"
24675   IF Block$="N" THEN GOTO Nobloin
24680  Nextb: PRINT "  B";Ib,
24685   Sgf=0                          ! Stopgap flag.
24690   Id=0
24695  Nextid: Id=Id+1
24700   GOSUB Subbloin
24705   IF Last=1 THEN Dinend          ! End of whole block
24710   IF Id<>Nd THEN GOTO Nextid     ! Next data
24715   BEEP 2000,1                    ! End of each block
24720   INPUT "Correction (Data #, NO=0) ?",Id
24725   IF Id=0 THEN
24730     GOTO 24755
24735   ELSE
24740     GOSUB Subbloin
24745     GOTO 24720
24750   END IF
24755   PRINT
24760   Ib=Ib+1
24765   GOTO Nextb
24770  Dinend: Nb=Ib-1
24775   Last=0
24780   IF Nr<>0 THEN REDIM Rdata(Nb,Nr)
24785   IF Ni<>0 THEN REDIM Idata(Nb,Ni)
24790   IF Nc<>0 THEN REDIM Cdata$(Nb,Nc)
24795   PRINT
24800   PRINT "Correction"
24805   INPUT "More correction (Block#,Data#,NO=0,0) ?",Ib,Id
24810   IF Ib<>0 THEN
24815     PRINT "  B";Ib,
24820     GOSUB Subbloin
24825     GOTO 24805
24830   END IF
24835   GOTO Out
24840  Nobloin:! ------------------Input unblocked data----------------
24845   I=1
24850   GOSUB Subnobloin
24855   IF Last=1 THEN 24870
24860   I=I+1
24865   GOTO 24850
24870   Nd=I-1
24875   Nr=0
24880   Ni=0
24885   Nc=0
24890   SELECT Item$(1,2)
24895   CASE "R"
24900     REDIM Rdata(1,Nd)
24905     Nr=Nd
24910   CASE "I"
24915     REDIM Idata(1,Nd)
24920     Ni=Nd
24925   CASE "C"
24930     REDIM Cdata$(1,Nd)
24935     Nc=Nd
24940   END SELECT
24945   PRINT
24950   Last=0
24955   PRINT "Correction"
24960   INPUT "Correction (Data#,NO=0) ?",I
24965   IF I<>0 THEN
24970     GOSUB Subnobloin
24975     GOTO 24960
24980   END IF
24985   ! ------------------Print-------------------------------
24990  Out: PRINT
24995   PRINT
25000   PRINT "Numbers of data"
25005   PRINT "Np","Nb","Nd","Nr","Ni","Nc"
25010   PRINT Np,Nb,Nd,Nr,Ni,Nc
25015   PRINT
25020   PRINT "Final data"
25025   IF Block$="N" THEN
25030     FOR Id=1 TO Nd
25035       PRINT Id;":";
25040       SELECT Item$(1,2)
25045       CASE "R"
25050         PRINT Rdata(1,Id),
25055       CASE "I"
25060         PRINT Idata(1,Id),
25065       CASE "C"
25070         PRINT Cdata$(1,Id),
25075       END SELECT
25080     NEXT Id
25085     PRINT
25090   ELSE
25095     FOR Ib=1 TO Nb
25100       PRINT "  B";Ib,
25105       FOR Id=1 TO Nd
25110         PRINT Id;":";
25115         SELECT Item$(Id,2)
25120         CASE "R"
25125           PRINT Rdata(Ib,Itable(Id)),
25130         CASE "I"
25135           PRINT Idata(Ib,Itable(Id)),
25140         CASE "C"
25145           PRINT Cdata$(Ib,Itable(Id)),
25150         END SELECT
25155       NEXT Id
25160       PRINT
25165     NEXT Ib
25170   END IF
25175   ! --------------Write & read file-----------------------
25180   Fileio("WR",Fn$,Comment$,Nb,Np,Nr,Ni,Nc,Paraname$(*),Paraval(*),Item$(*),Rdata(*),Idata(*),Cdata$(*))
25185   PRINTER IS CRT
25190   SUBEXIT
25195   ! --------------------SUBBloin--------------------------
25200  Subbloin:!
25205   Dispseq$="))"
25210   IF Seq(Id,1)=1 THEN Dispseq$=","&VAL$(Seq(Id,2))&":Seq.data end))"
25215   SELECT Item$(Id,2)
25220   CASE "R"
25225     DISP "Block ";Ib;"  Data ";Id;" (";Item$(Id,1);" ,R (-1234:whole end";Dispseq$,
25230     Ir=Itable(Id)
25235     INPUT Rdata(Ib,Ir)
25240     IF Seq(Id,1)=1 AND Rdata(Ib,Ir)=Seq(Id,2) THEN GOSUB Endp
25245     IF Rdata(Ib,Ir)=-1234 THEN Last=1
25250     IF Last<>1 AND Seq(Id,2)<>Rdata(Ib,Ir) THEN PRINT Id;":";Rdata(Ib,Ir),
25255   CASE "I"
25260     DISP "Block ";Ib;"  Data ";Id;" (";Item$(Id,1);" ,I, -1234 to block end)",
25265     Ii=Itable(Id)
25270     INPUT Idata(Ib,Ii)
25275     IF Idata(Ib,Ii)=Seq(Id,2) THEN GOSUB Endp
25280     IF Idata(Ib,Ii)=-1234 THEN Last=1
25285     IF Last<>1 AND Seq(Id,2)<>Idata(Ib,Ii) THEN PRINT Id;":";Idata(Ib,Ii),
25290   CASE "C"
25295     DISP "Block ";Ib;"  Data ";Id;" (";Item$(Id,1);" ,C, -1234 to block end)",
25300     Ic=Itable(Id)
25305     LINPUT Cdata$(Ib,Ic)
25310     IF Cdata$(Ib,Ic)=VAL$(Seq(Id,2)) THEN GOSUB Endp
25315     IF Cdata$(Ib,Ic)="-1234" THEN Last=1
25320     IF Last<>1 AND VAL$(Seq(Id,2))<>Cdata$(Ib,Ic) THEN PRINT Id;":";Cdata$(Ib,Ic),
25325   END SELECT
25330   RETURN
25335  Subnobloin: ! ------------------SUBNobloin--------------------------
25340   DISP "Data ";I;" (-1234 to end)",
25345   SELECT Item$(1,2)
25350   CASE "R"
25355     INPUT Rdata(1,I)
25360     IF Rdata(1,I)=-1234 THEN Last=1
25365     IF Last<>1 THEN PRINT I;":";Rdata(1,I),
25370   CASE "I"
25375     INPUT Idata(1,I)
25380     IF Idata(1,I)=-1234 THEN Last=1
25385     IF Last<>1 THEN PRINT I;":";Idata(1,I),
25390   CASE "C"
25395     INPUT Cdata$(1,I)
25400     IF Cdata$(1,I)="-1234" THEN Last=1
25405   END SELECT
25410   RETURN
25415  Inc: ! ------------------------------------------------------
25420   SELECT Item$(Id,2)
25425   CASE "R"
25430     Ir=Ir+1
25435     Itable(Id)=Ir
25440   CASE "I"
25445     Ii=Ii+1
25450     Itable(Id)=Ii
25455   CASE "C"
25460     Ic=Ic+1
25465     Itable(Id)=Ic
25470   END SELECT
25475   Id=Id+1
25480   RETURN
25485  Endp:     ! ------------------------------------------------------
25490   SELECT Item$(Id,2)
25495   CASE "R"
25500     Rdata(Ib,Ir)=Seq(Id,2)
25505     PRINT Id;":";Rdata(Ib,Ir),
25510     Id=Id+1
25515     Ir=Ir+1
25520     IF Id<=Nd AND Ir<=Nr THEN
25525       IF Seq(Id,1)=1 THEN
25530         GOTO 25500
25535       ELSE
25540         Id=Id-1
25545         Ir=Ir-1
25550         RETURN
25555       END IF
25560     ELSE
25565       Id=Id-1
25570       Ir=Ir-1
25575       RETURN
25580     END IF
25585   CASE "I"
25590     Idata(Ib,Ii)=Seq(Id,2)
25595     PRINT Id;":";Idata(Ib,Ii),
25600     Ii=Ii+1
25605     Id=Id+1
25610     IF Id<=Nd AND Ii<=Ni THEN
25615       IF Seq(Id,1)=1 THEN
25620         GOTO 25590
25625       ELSE
25630         Id=Id-1
25635         Ii=Ii-1
25640         RETURN
25645       END IF
25650     ELSE
25655       Id=Id-1
25660       Ii=Ii-1
25665       RETURN
25670     END IF
25675   CASE "C"
25680     Cdata$(Ib,Ic)=VAL$(Seq(Id,2))
25685     PRINT Id;":";Cdata$(Ib,Ic),
25690     Ic=Ic+1
25695     Id=Id+1
25700     IF Id<=Nd AND Ic<=Nc THEN
25705       IF Seq(Id,1)=1 THEN
25710         GOTO 25680
25715       ELSE
25720         Id=Id-1
25725         Ic=Ic-1
25730         RETURN
25735       END IF
25740     ELSE
25745       Id=Id-1
25750       Ic=Ic-1
25755       RETURN
25760     END IF
25765   END SELECT
25770   RETURN
25775 SUBEND
25780 ! ####################################################
25785 SUB Fileio(Wr$,Filename$,Comment$,INTEGER Nb,Np,Nr,Ni,Nc,Paraname$(*),Paraval(*),Item$(*),Rdata(*),INTEGER Idata(*),Cdata$(*))
25790 ! Wr$  R:Read only, W:Write only, RW:Write & Read
25795 ! Data are written or read in order of real,int, and cha.
25800 ! Exam. of DIM in calling pro.(COPY & CHANGE). These are redimed.
25805 !    DIM Filename$[10],Com$[100],Para$(20)[20],Para(20),Item$(100,2)[10],Rd(10,100),Cd$(1,1)[10]
25810 !    INTEGER Id(1,1),Nk,Np,Nr,Ni,Nc
25815 ! ####################################################
25820   OPTION BASE 1
25825   Filename(Filename$)
25830   IF Wr$="R" THEN GOTO Rf
25835   Trl=(4+10)+(4+100)+(2*5)+((4+20)*Np)+(8*Np)+(2*(4+10)*SIZE(Item$,1))+((8*Nr)+(2*Ni)+((4+10)*Nc))*Nb+1000
25840   Rl=256
25845   Rn=INT(Trl/Rl)+1
25850   PRINT
25855   PRINT "Memory needed",Rn*Rl/1000;"[kByte]"
25860 ! ------------------Write to file-----------------------
25865   CREATE BDAT Filename$,Rn,Rl
25870   DISP "Writing"
25875   ASSIGN @File TO Filename$
25880   OUTPUT @File;Filename$,Comment$
25885   OUTPUT @File;Nb,Np,Nr,Ni,Nc
25890   IF Np<>0 THEN OUTPUT @File;Paraname$(*),Paraval(*)
25895   OUTPUT @File;Item$(*)
25900 ! ------------------Unblocked data----------------------
25905   IF Nb=1 THEN
25910     SELECT Item$(1,2)
25915     CASE "R"
25920       OUTPUT @File;Rdata(*),END
25925     CASE "I"
25930       OUTPUT @File;Idata(*),END
25935     CASE "C"
25940       OUTPUT @File;Cdata$(*),END
25945     END SELECT
25950 ! ------------------Blockd data-------------------------
25955   ELSE
25960     FOR Ib=1 TO Nb
25965       IF Nr<>0 THEN
25970         FOR Ir=1 TO Nr
25975           OUTPUT @File;Rdata(Ib,Ir),END
25980         NEXT Ir
25985       END IF
25990       IF Ni<>0 THEN
25995         FOR Ii=1 TO Ni
26000           OUTPUT @File;Idata(Ib,Ii),END
26005         NEXT Ii
26010       END IF
26015       IF Nc<>0 THEN
26020         FOR Ic=1 TO Nc
26025           OUTPUT @File;Cdata$(Ib,Ic),END
26030         NEXT Ic
26035       END IF
26040     NEXT Ib
26045   END IF
26050   BEEP
26055   PRINT
26060   PRINT "Written to file:";Filename$
26065   IF Wr$="W" THEN SUBEXIT
26070  Rf: ! ------------------Read from file----------------------
26075   ASSIGN @File TO Filename$
26080   DISP "Reading"
26085   ENTER @File;Filenamep$,Comment$
26090   IF Filenamep$<>Filename$ THEN
26095     BEEP
26100     DISP "Filename is not matched. Make sure & <CONT>."
26105     PAUSE
26110     GOTO 25825
26115   END IF
26120   ENTER @File;Nb,Np,Nr,Ni,Nc
26125   Nd=Nr+Ni+Nc
26130   IF Np<>0 THEN
26135     REDIM Paraname$(Np),Paraval(Np)
26140     ENTER @File;Paraname$(*),Paraval(*)
26145   END IF
26150   ! ------------------------------------------------------
26155   SELECT Nb
26160   CASE 1
26165     REDIM Item$(1,2)
26170     ENTER @File;Item$(*)
26175   CASE ELSE
26180     REDIM Item$(Nd,2)
26185     ENTER @File;Item$(*)
26190     ALLOCATE Itable(Nd)
26195     Ir=0
26200     Ii=0
26205     Ic=0
26210     FOR Id=1 TO Nd
26215       SELECT Item$(Id,2)
26220       CASE "R"
26225         Ir=Ir+1
26230         Itable(Id)=Ir
26235       CASE "I"
26240         Ii=Ii+1
26245         Itable(Id)=Ii
26250       CASE "C"
26255         Ic=Ic+1
26260         Itable(Id)=Ic
26265       END SELECT
26270     NEXT Id
26275   END SELECT
26280   ! ------------------------------------------------------
26285   IF Nr<>0 THEN REDIM Rdata(Nb,Nr)
26290   IF Ni<>0 THEN REDIM Idata(Nb,Ni)
26295   IF Nc<>0 THEN REDIM Cdata$(Nb,Nc)
26300   FOR Ib=1 TO Nb
26305     IF Nb=1 AND Ib<>1 THEN GOTO 26395
26310     IF Nr<>0 THEN
26315       FOR Ir=1 TO Nr
26320         ENTER @File;Rdata(Ib,Ir)
26325       NEXT Ir
26330     END IF
26335     IF Ni<>0 THEN
26340       FOR Ii=1 TO Ni
26345         ENTER @File;Idata(Ib,Ii)
26350       NEXT Ii
26355     END IF
26360     IF Nc<>0 THEN
26365       FOR Ic=1 TO Nc
26370         ENTER @File;Cdata$(Ib,Ic)
26375       NEXT Ic
26380     END IF
26385   NEXT Ib
26390      ! ------------------Print-------------------------------
26395   INPUT "Printer(C/P/NO) ?",Prt$
26400   IF Prt$="NO" THEN SUBEXIT
26405   IF Prt$="P" THEN PRINTER IS 701
26410   PRINT
26415   PRINT "**************************************************"
26420   PRINT "Read from file:";Filename$
26425   PRINT
26430   PRINT "Comment"
26435   PRINT TAB(2),Comment$
26440   IF Np<>0 THEN
26445     PRINT
26450     PRINT "Parameters"
26455     FOR Ip=1 TO Np
26460       PRINT TAB(2),Paraname$(Ip),Paraval(Ip)
26465     NEXT Ip
26470   END IF
26475   ! ------------------------------------------------------
26480   PRINT
26485   SELECT Nb
26490   CASE 1
26495     PRINT "Data name & Type",Item$(1,1),Item$(1,2)
26500     PRINT "  ";Item$(1,1),Item$(1,2)
26505   CASE ELSE
26510     PRINT "Data name & Type"
26515     PRINT "No.","Name","Type","No. in type"
26520     FOR Id=1 TO Nd
26525       PRINT Id,Item$(Id,1),Item$(Id,2),Itable(Id)
26530     NEXT Id
26535   END SELECT
26540   PRINT
26545   PRINT "Numbers of data"
26550   PRINT "Np","Nb","Nd","Nr","Ni","Nc"
26555   PRINT Np,Nb,Nd,Nr,Ni,Nc
26560   ! ------------------------------------------------------
26565   PRINT
26570   PRINT "Data"
26575   IF Nb=1 THEN
26580     FOR Id=1 TO Nd
26585       PRINT Id;":";
26590       SELECT Item$(1,2)
26595       CASE "R"
26600         PRINT Rdata(1,Id),
26605       CASE "I"
26610         PRINT Idata(1,Id),
26615       CASE "C"
26620         PRINT Cdata$(1,Id),
26625       END SELECT
26630     NEXT Id
26635     PRINT
26640   ELSE
26645     FOR Ib=1 TO Nb
26650       PRINT "  B";Ib,
26655       FOR Id=1 TO Nd
26660         PRINT Id;":";
26665         SELECT Item$(Id,2)
26670         CASE "R"
26675           PRINT Rdata(Ib,Itable(Id)),
26680         CASE "I"
26685           PRINT Idata(Ib,Itable(Id)),
26690         CASE "C"
26695           PRINT Cdata$(Ib,Itable(Id)),
26700         END SELECT
26705       NEXT Id
26710       PRINT
26715     NEXT Ib
26720   END IF
26725   IF Prt$="P" THEN PRINTER IS CRT
26730 SUBEND
26735 SUB Tsavemeas
26740 ! NOF averaging of measured TS functions.
26745 ! Format in Filein
26750 !  Blocking mode. If one block, block 2 must be dummy.
26755 !  Para  1:Stopgap=200, 2:Ang.step, 3:Max.ang.
26760 !  Idata 1:#, 2:Year, 3:Month, 4:Day, 5:Time, 6:Freq.
26765 !  Cdata 7:Species, 8:Comment
26770 !  Rdata 9:Length, 10:Height, 11:Breadth, 12:Weight, 13 to end:-TS
26775 ! ####################################################
26780   OPTION BASE 1
26785   DIM Para(3),Para$(3)[10],Item$(170,2)[10],Cd$(5,2)[10],Rd(5,150),Fn$[10],Com$[100]
26790   INTEGER Id(5,6),Nb,Np,Nr,Ni,Nc
26795   Fileio("R",Fn$,Com$,Nb,Np,Nr,Ni,Nc,Para$(*),Para(*),Item$(*),Rd(*),Id(*),Cd$(*))
26800   ALLOCATE Ts(Nr-4)
26805   FOR I=1 TO Nr-4
26810     Ts(I)=-Rd(1,I+4)
26815   NEXT I
26820   PRINT
26825   INPUT "Theta: Mean, Sigma ?",Meant,Sigmat
26830   PRINT "Theta: Mean, Sigma",Meant,Sigmat
26835   Tsa=FNTsnof(Ts(*),Para(2),Meant,Sigmat)
26840   PRINT "Ave. TS",Tsa;"dB"
26845   INPUT "Other orient.(O), TS func.(T) or exit(E) ?",Oce$
26850   SELECT Oce$
26855   CASE "O"
26860     GOTO 26825
26865   CASE "T"
26870     GOTO 26795
26875   CASE "E"
26880     SUBEXIT
26885   END SELECT
26890 SUBEND
26895 ! ####################################################
26900 SUB Polax(Tmin,Tmax,Torg,Rmin,Rmax,Axlas,Rlb$,Ttl$,Auto)
26905 ! Axes, scale, and label for polar graph.
26910 ! T's are angles refered to +x axis. Torg is origin of labeled axis.
26915 ! If Axlas=0 then only write R_axis value when Axl<>360.( See line Axl )
26920 ! Origin is (0,0) & radial max. scale is Rsmax-Rsmin.
26925 ! If Rlb$=" " or Plot$="705" then labels can be keyed in.
26930 ! Rmin & max is changed to Rsmin & max.
26935 ! ####################################################
26940   COM /Gdu/Xgmax,Ygmax,Upg,Upg_      ! Upg ->       Upg_:dummy
26945   COM /Margin/Xlm,Xrm,Ybm,Ytm        !     <-
26950   COM /Range/Xsmin,Xsmax,Xsr,Ysmin,Ysmax,Ysr,Csl   !  ->
26955   COM /Auto/Pltno,Tg$,Pen$,Tgs,Tts,Rsmin,Rsmax,Rgs,Rts,Axl,Line_pen$,Posi$,Mf$,Filename$   ! <--
26960   COM /Pltr/Plot$[8],Sp$,Lastpen             !     <-
26965   DIM Ax$[4],Axl$[8]
26970   DEG
26975   Rln=35                     ! Maximum character number.
26980   Tln=45
26985   ALLOCATE Rl$[Rln+50],Tytle$[Tln+50] ! 50 is margin for input
26990   CALL Graphinit(1,Auto)
26995   ! ------------------Input-------------------------------
27000   Rl$=Rlb$
27005   Tytle$=Ttl$
27010   IF Rlb$=" " OR Plot$="705" THEN
27015     INPUT "Keyed label ?(Y/N)",Kl$
27020     IF Kl$="Y" THEN
27025       LINPUT "R label ?",Rl$
27030       LINPUT "Tytle ?",Tytle$
27035     END IF
27040   END IF
27045   CALL Label_length(Rl$,Rln)
27050   CALL Label_length(Tytle$,Tln)
27055   IF Auto=0 THEN
27060     INPUT "Tic(T) or grid(G) ?",Tg$
27065     INPUT "Do you want two pens ? (Y/N)",Pen$
27070     DISP "Rmin=";Rmin;", Rmax=";Rmax;"     ";! Not auto, so you can
27075     INPUT "R_scale_min, _max ?",Rsmin,Rsmax  ! specify plotting range.
27080     Rsr=Rsmax-Rsmin
27085     DISP "Rsmin=";Rsmin;", Rsmax=";Rsmax;", R_range=";Rsr;"     ";
27090     INPUT "R_grid_&_label_space, R_tic_space ?",Rgs,Rts
27095     Tsr=Tmax-Tmin
27100     DISP "Tmin=";Tmin;", Tmax=";Tmax;", Tsr=";Tsr;"     ";
27105     INPUT "Angle_grid_, and Tic_space ? (No:0,0)",Tgs,Tts
27110  Axl: INPUT "To which axes do you want to label ? (Exam.90,No=360)",Axl
27115          ! If Axlas=0 then only scale values are labeled.
27120   ELSE
27125     Rsr=Rsmax-Rsmin
27130     Tsr=Tmax-Tmin
27135   END IF
27140   ! ------------------Scaling-----------------------------
27145   Rmin=Rsmin
27150   Rmax=Rsmax
27155   Xsmin=1
27160   Xsmax=-1
27165   Ysmin=1
27170   Ysmax=-1
27175   FOR T=Tmin TO Tmax STEP (Tmax-Tmin)/20  ! Find range.
27180     St=SIN(T)
27185     Ct=COS(T)
27190     IF St>Ysmax THEN Ysmax=St
27195     IF St<Ysmin THEN Ysmin=St
27200     IF Ct>Xsmax THEN Xsmax=Ct
27205     IF Ct<Xsmin THEN Xsmin=Ct
27210   NEXT T
27215   IF Xsmax<0 THEN Xsmax=0       ! Origin must be in range.
27220   IF Xsmin>0 THEN Xsmin=0
27225   IF Ysmax<0 THEN Ysmax=0
27230   IF Ysmin>0 THEN Ysmin=0
27235   Marg=.2
27240   Xsmin=Rsr*(Xsmin-Marg)
27245   Xsmax=Rsr*(Xsmax+Marg)
27250   Ysmin=Rsr*(Ysmin-(Axlas=1)*Marg-Marg)   ! Tytle margin.
27255   Ysmax=Rsr*(Ysmax+Marg)
27260   Xsr=Xsmax-Xsmin
27265   Ysr=Ysmax-Ysmin
27270   SHOW Xsmin,Xsmax,Ysmin,Ysmax
27275   ! --------------------Factors---------------------------
27280   Xgw=Xgmax*(1-(Xlm+Xrm)/100)    ! X gdu width.
27285   Ygw=Ygmax*(1-(Ybm+Ytm)/100)
27290   Xbyyg=Xgw/Ygw
27295   Xbyys=Xsr/Ysr
27300   IF Xbyyg>Xbyys THEN
27305     Ymg=0
27310     Xmg=(Ysr*Xbyyg-Xsr)/2
27315   ELSE
27320     Xmg=0
27325     Ymg=(Xsr/Xbyyg-Ysr)/2
27330   END IF
27335   Xsmin=Xsmin-Xmg
27340   Xsmax=Xsmax+Xmg
27345   Xsr=Xsmax-Xsmin
27350   Ysmin=Ysmin-Ymg
27355   Ysmax=Ysmax+Ymg
27360   Ysr=Ysmax-Ysmin
27365   Upg=Xsr/Xgw
27370   Upg_=Ysr/Ygw
27375   Csv=Ysr/29/Upg                  ! Udu value is used to make size
27380   Csl=Ysr/25/Upg                  ! proportionate to range.
27385   Cst=Ysr/22/Upg
27390   Ticlen=Rsr/70
27395   Tic2=Ticlen*2
27400   ! --------------------Axes and R_tic--------------------
27405   Ax$=""                        ! Find Axes.
27410   FOR T=0 TO 270 STEP 90
27415     IF T>=Tmin AND T<=Tmax THEN
27420       Ax$=Ax$&"1"               ! Exam. Ax$="1010" means 90 & 270.
27425     ELSE
27430       Ax$=Ax$&"0"
27435     END IF
27440   NEXT T
27445   FOR I=0 TO 3                  ! Draw axes.
27450     IF Ax$[I+1;1]="0" THEN Nexti
27455     MOVE 0,0
27460     PIVOT I*90
27465     Lmax=1
27470     IF I=0 THEN Lmax=5
27475     IF I=1 THEN Lmax=2
27480     FOR L=1 TO Lmax
27485       DRAW Rsr,0
27490       MOVE 0,0
27495     NEXT L
27500     FOR Decade=Rsmin TO Rsmax-Rgs STEP Rgs   ! Draw tic.
27505       FOR Units=1-(Decade=Rsmin) TO Rgs/Rts
27510         R=Decade-Rsmin+Units*Rts
27515         Tl=Ticlen*(1+(Units=Rgs/Rts))
27520         MOVE R,-Tl
27525         DRAW R,Tl
27530       NEXT Units
27535     NEXT Decade
27540     PIVOT 0
27545  Nexti: NEXT I
27550   ! ------------------Arc(R_grid)-------------------------
27555   IF Tgs=0 THEN Rav
27560   Arcdeg(0,0,11,0,Rsr,Tmin,Tsr)  ! Outer-most arc.
27565   IF Tg$="G" THEN
27570     FOR R=Rsmin+Rgs TO Rsmax STEP Rgs  ! Shifted to avoid overlap.
27575       Arcdeg(0,0,11,0,R-Rsmin,Tmin,Tsr)
27580     NEXT R
27585   END IF
27590   ! ------------------Angle tic---------------------------
27595   FOR Decade=Tmin TO Tmax STEP Tgs  ! Draw angle tic.
27600     Dg=((Decade=Tmin AND (Tmin MOD 90)<>0) OR (Decade=Tmax AND (Tmax MOD 90)<>0))                                       ! End lines
27605     IF Dg=1 OR ((Decade MOD 90<>0) AND Tg$="G") THEN
27610       MOVE 0,0
27615       PIVOT Decade
27620       DRAW Rsr,0
27625     END IF
27630     FOR Units=1 TO Tgs/Tts
27635       MOVE 0,0
27640       T=Decade+Units*Tts
27645       PIVOT T
27650       IF (Tg$="T" AND T<Tmax) OR (Tg$="G" AND Units<>Tgs/Tts) THEN
27655         MOVE Rsr-Ticlen*(1+(Units=Tgs/Tts)),0
27660         DRAW Rsr,0
27665       END IF
27670     NEXT Units
27675   NEXT Decade
27680   PIVOT 0
27685  Rav:! --------------------R axis value----------------------
27690   IF Pen$="Y" THEN PEN 2
27695   CSIZE Csv
27700   SELECT Axl
27705   CASE 0,180
27710     FOR R=Rsmin TO Rsmax STEP Rgs
27715       LORG 6
27720       Rr=R-Rsmin
27725       IF Axl=180 THEN Rr=-Rr
27730       MOVE Rr,-Tic2
27735       LABEL USING "#,K";R
27740     NEXT R
27745   CASE 90,270
27750     LORG 2
27755     Maxlen=1                        ! For R label.
27760     FOR R=Rsmin TO Rsmax STEP Rgs
27765       Length=LEN(VAL$(R))
27770       IF Length>Maxlen THEN Maxlen=Length
27775     NEXT R
27780     Csvx=Csv*Upg*9/15
27785     Xv=-Csvx*Maxlen-Ticlen*3
27790     FOR R=Rsmin TO Rsmax STEP Rgs
27795       Rr=R-Rsmin
27800       IF Axl=270 THEN Rr=-Rr
27805       MOVE Xv+(Maxlen-1-LEN(VAL$(R)))*Csvx,Rr
27810       LABEL USING "#,X,K";R
27815     NEXT R
27820   END SELECT
27825   ! --------------------Angle value----------------------
27830   IF Tgs=0 THEN Rl
27835   IF Axlas=0 THEN SUBEXIT
27840   LORG 5
27845   CSIZE Csv
27850   Rr=Rsr+Csv*Upg
27855   FOR T=Tmin TO Tmax-(Tmax=360)*Tgs STEP Tgs
27860     MOVE Rr*COS(T),Rr*SIN(T)
27865     LDIR T-90
27870     LABEL USING "#,K";T-Torg
27875   NEXT T
27880  Rl:!--------------------R label---------------------------
27885   CSIZE Csl
27890   Rly=-Tic2-Csv*Upg*1.2
27895   Rr=Rsr/2-Csl*Upg*LEN(Rl$)*9/15/2
27900   SELECT Axl
27905   CASE 0,180
27910     LORG 6
27915     LDIR 0
27920     IF Axl=0 THEN MOVE Rr,Rly
27925     IF Axl=180 THEN MOVE -Rsr+Rr,Rly
27930   CASE 90,270
27935     Rxlp=Xv
27940     LORG 1
27945     LDIR 90
27950     IF Axl=90 THEN MOVE Rxlp,Rr
27955     IF Axl=270 THEN MOVE Rxlp,Rsr-Rr
27960   CASE 360                           ! Not write. See line Axl.
27965     GOTO Tytle
27970   END SELECT
27975   Label(Rl$)
27980  Tytle:! ------------------Tytle-------------------------------
27985   LORG 4
27990   LDIR 0
27995   CSIZE Cst
28000   Txp=(Xsmax+Xsmin)/2-Cst*Upg*LEN(Tytle$)*9/15/2
28005   Typ=Ysmin+Cst*Upg
28010   MOVE Txp,Typ
28015   Label(Tytle$)
28020 SUBEND
28025   ! ####################################################
28030 DEF FNSjbes(Nn,Xx)
28035 ! ?/2118
28040 ! Spherical Bessel function of first kind.
28045 ! See "Value Analysis and FORTRAN" by Amamiya and Taguchi,
28050 !     Maruzen Co.,Tokyo,1971.
28055 ! #######################(SF-5)#######################
28060   Dx0=.0000003                ! use to aboid X=integer*PI
28065   N=Nn
28070   X=Xx
28075   IF N>30000 THEN GOSUB Notaccu
28080   SELECT X
28085   CASE <0
28090     GOSUB Invalid
28095   CASE 0 TO 7.E-4
28100     GOSUB Small
28105   CASE 7.E-4 TO 3.E+4
28110     GOSUB Large
28115   CASE >3.E+4
28120     GOSUB Notaccu
28125   END SELECT
28130 ! ------------------------------------------------------
28135  Invalid: PRINT "Argument of SJBES is invalid.  N=";N;"  X=";X
28140   RETURN 0
28145  Notaccu: PRINT "Value of SJBES is not accurate.  N=";N;"  X=";X
28150   RETURN 0
28155 ! ------------------------------------------------------
28160  Large: SELECT X
28165   CASE <.2
28170     Y=X*X
28175     W=1-Y*(1-.05*Y)/6          ! Approximation of sinX/X.
28180   CASE >=.2
28185     M=INT(X/PI)
28190     Dx=X-M*PI
28195     SELECT Dx
28200     CASE <Dx0
28205       X=M*PI+Dx0              ! Aboid integer*PI
28210     CASE >PI-Dx0
28215       X=(M+1)*PI-Dx0
28220     CASE ELSE
28225       X=X
28230     END SELECT
28235     W=FNSinx(X)/X            ! Can't use SIN for deep call.
28240   END SELECT
28245   SELECT N
28250   CASE <0
28255     GOSUB Invalid
28260   CASE 0
28265     RETURN W
28270   END SELECT
28275   SELECT X
28280   CASE >=100
28285     L=.02*X+18
28290   CASE 10 TO 100
28295     L=.1*X+10
28300   CASE 1 TO 10
28305     L=.5*X+5
28310   CASE <1
28315     L=5
28320   END SELECT
28325   Nm=MAX(N,INT(X))+INT(L)
28330   Z=1/X
28335   T3=0
28340   T2=1.E-35
28345   FOR I=1 TO Nm
28350     K=Nm-I
28355     T1=(K+K+3)*Z*T2-T3
28360     IF ABS(N-K)<.01 THEN Sj=T1
28365     IF ABS(T1)>=1.E+25 THEN
28370       T1=T1*1.E-25
28375       T2=T2*1.E-25
28380       Sj=Sj*1.E-25
28385     END IF
28390     T3=T2
28395     T2=T1
28400   NEXT I
28405   RETURN W/T1*Sj
28410 ! ------------------------------------------------------
28415  Small: W=1
28420   SELECT N
28425   CASE <0
28430     GOSUB Invalid
28435   CASE 0
28440     RETURN W
28445   CASE 0 TO 10
28450     T1=3
28455     T2=1
28460     FOR I=1 TO N
28465       T3=T2*X/T1
28470       T1=T1+2
28475       T2=T3
28480     NEXT I
28485     RETURN T3
28490   CASE >10
28495     RETURN 0
28500   END SELECT
28505 FNEND
28510   ! ####################################################
28515 DEF FNSinx(Y)
28520 ! ?/2118
28525 ! sin x to be used when HP SIN could not be used.
28530 ! ####################################################
28535   Accu=1.E-4
28540   Loopn=500
28545   X=Y
28550 ! IF X=0 THEN X=1.E-30        ! is necessary ? ######
28555   B=X
28560   S=X
28565   S1=0
28570   X2=X*X
28575   FOR D=2 TO Loopn STEP 2
28580     B=-B*X2/D/(D+1)
28585     S=S+B
28590     IF ABS((S-S1)/S)<Accu THEN RETURN S
28595     S1=S
28600   NEXT D
28605   BEEP
28610   PRINT "Accu in FNSin not satisfied."
28615   RETURN S
28620 FNEND
28625 ! ####################################################
28630 SUB Pltfile(Autor)
28635 ! */1218/2121
28640 ! Read, print, trasfer, & plot filed data made by Plot.
28645 !   Autor: 0 for manual, 1 for auto.
28650 ! ####################################################
28655   DIM Filename$[10],Rfilename$[10]
28660   COM /Pltd/Lnum,Dnum,Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,Auto,Autled$,Linlbl$(*),Manuled
28665   COM /Seqrpf/Basefname$[8],Fnmin,Fnmax,Op$[8],Pr$[8],Py$[8],Pl$[8],Fnum,Endflg
28670  Start:!----------------Input from file---------------------
28675   IF Autor=0 THEN
28680     Filename(Filename$)
28685   ELSE
28690     IF Fnum=Fnmax-Fnmin+1 THEN
28695       Endflg=1
28700       SUBEXIT
28705     END IF
28710     F2$=VAL$(Fnmin+Fnum)
28715     IF LEN(F2$)=1 THEN F2$="0"&F2$
28720     Filename$=Basefname$&F2$
28725     Fnum=Fnum+1
28730   END IF
28735   ASSIGN @File TO Filename$
28740   ON ERROR GOTO Check
28745   ENTER @File;Rfilename$
28750   OFF ERROR
28755   IF Filename$<>Rfilename$ THEN
28760     BEEP
28765     DISP "File name is not matched. Examine disc or file name & <CONT>."
28770     PAUSE
28775     Fnum=Fnum-1
28780     GOTO 28675
28785   END IF
28790   ENTER @File;Lnum,Dnum
28795   REDIM Y(Lnum-1,Dnum-1),Linlbl$(Lnum-1)
28800   ENTER @File;Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,Autled$,Linlbl$(*),Manuled
28805   ASSIGN @File TO *
28810   ! ------------------Print data---------------------------
28815   IF Autor=0 THEN INPUT "Print data(Y/N) ?",Op$
28820   IF Op$="N" THEN GOTO 29140
28825   IF Autor=0 THEN INPUT "Printer (C/P) ?",Pr$
28830   IF Pr$="P" THEN PRINTER IS 701
28835   PRINT
28840   PRINT "DATA READ FROM PLOTFILE:";Filename$
28845   PRINT "***********************************"
28850   PRINT
28855   PRINT "(1)Tytle",TAB(25),Ttl$
28860   PRINT "(2)Auto_legend",TAB(25),Autled$
28865   PRINT "(3)Y_Label",TAB(25),Ylb$
28870   PRINT "(4)X_Label",TAB(25),Xlb$
28875   PRINT "(5-6)X_min, _max, _inc",TAB(25),Xmin,Xmax,Xinc
28880   PRINT "Line & Data number",TAB(25),Lnum,Dnum
28885   PRINT "Line label (LL)";
28890   FOR L=0 TO Lnum-1
28895     PRINT "(";L+7;")",TAB(25),L,Linlbl$(L)
28900   NEXT L
28905   PRINT "(";Lnum+7;")Graph kind",TAB(25),G$
28910   PRINT "(";Lnum+8;")Manual legend",TAB(25),Manuled
28915   IF Autor=0 THEN
28920     INPUT "Change parameter (Y/N) ?",Cp$
28925     IF Cp$="Y" THEN
28930       INPUT "Which para. (by number) ?",Pnum
28935       SELECT Pnum
28940       CASE 6
28945         INPUT "X_max ?",Xmax
28950       CASE ELSE
28955         PRINT "Not cmpleted."
28960       END SELECT
28965       INPUT "More change (Y/N) ?",Mc$
28970       IF Mc$="Y" THEN GOTO 28930
28975     END IF
28980     INPUT "Print Y data (Y/N) ?",Py$
28985     IF Py$="N" THEN GOTO 29140
28990   END IF
28995   PRINT
29000   PRINT "Data#          X";
29005   FOR L=0 TO Lnum-1
29010     Form$="#,7X,""LL "",D"
29015     PRINT USING Form$;L
29020   NEXT L
29025   FOR D=0 TO Dnum-1
29030     PRINT
29035       SELECT G$
29040       CASE "LIN"
29045         X=Xmin+D*Xinc
29050       CASE "POL"
29055         X=Xmin+D*Xinc
29060 !       IF X>Xmax THEN GOTO 17000
29065       CASE "LOG","SML"
29070         IF Xinc>0 THEN
29075           X=Xmin*Xinc^D
29080         ELSE
29085           X=Xmin-D*Xinc
29090         END IF
29095       END SELECT
29100     PRINT USING "#,2X,3D,2X,SD.DDE";D,X
29105     FOR L=0 TO Lnum-1
29110       PRINT USING "#,2X,SD.DDE";Y(L,D)
29115     NEXT L
29120   NEXT D
29125   PRINT
29130   IF Pr$="P" THEN PRINTER IS CRT
29135   ! ------------------Plot--------------------------------
29140   IF Autor=0 THEN INPUT "Plot (Y/N) ?",Pl$
29145   IF Pl$="N" THEN SUBEXIT
29150   Plot(Y(*),Xmin,Xmax,Xinc,G$,Xlb$,Ylb$,Ttl$,0,Autled$,Linlbl$(*),1)
29155   IF F2$=Fnmax$ THEN SUBEXIT
29160   SUBEXIT
29165  Check:!---------------------------------------------------
29170   DISP "Check file name."
29175   GOTO Start
29180   RETURN
29185 SUBEND
29190 ! ####################################################
29195 DEF FNGamma(M,R,C)
29200 !  See Flammer eq.(3.1.5).
29205 ! ####################################################
29210   IF R<0 THEN
29215     BEEP
29220     OUTPUT 1;"r>=0 must be satisfied."
29225   END IF
29230   Mr=M+R
29235   Mr2=2*Mr
29240   RETURN Mr*(Mr+1)+C^2/2*(1-(4*M^2-1)/(Mr2-1)/(Mr2+3))
29245 FNEND
29250 !*  29250 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29260 !*  29260 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29270 !*  29270 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29280 !*  29280 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29290 !*  29290 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29300 !*  29300 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29310 !*  29310 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29320 !*  29320 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29330 !*  29330 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29340 !*  29340 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29350 !*  29350 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29360 !*  29360 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29370 !*  29370 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29380 !*  29380 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
29390 !*  29390 CSUB Gdump_colored(From_ds,To_ds,OPTIONAL Rotate$,INTEGER Resolution,Background$,Algorithm$)
