XProfan

WindowTitle "Unvollständige Elliptische Integrale 1. und 2. Art {Legendre-F(phi,k), E(phi,k)}"
' Versuch (Early pre-alpha) - ohne jegliche Gewähr!
'{Q: F77-Bibliotheksmodul 1982, HP41C UPN-Programm 1983 und TrueBasic-Modul 1997
' 2015-07 nach XProfan 11.2a montiert - ACHTUNG: Ausschließlich als Demo!
' *** Möglicherweise bestehen Rechte Dritter! *** MfG Z. Ahnpasta
Windowstyle 24:Cls:font 2:set("decimals",17)
declare q!,k!,a!,d2r! , F!,Ellpci_e! : d2r!=pi()/180

Whileloop 0,900,1'6

    q!=d2r!* &Loop/10'30 '°
    k!=sin(q!)'Hinweis: Das ist NICHT die 'Modulwinkel'-Version der Legendre-Normalform,
    'k soll hier einfach nur den Definitionsbereich [0..1] durchlaufen!
    print "\n k = ";format$("00.0",&Loop);" °"
    Waitinput 2000' oder auskommentieren ...

    whileloop 0,90,15 : a!=d2r!* &Loop/10'° in Zehntelgrad-Schritten

        '*************************************************************************
        'Proc-Aufruf setzt 2 Ergebnisse:
        F!=Ellpci(a!,k!)
        'Die ausserhalb des Proc definierte Variable Ellpci_E!
        'bringt hier das 2. Ergebnis, E(phi,k) (Legendre-Form 2. Art) zurück.
        '*************************************************************************
        '*************************************************************************
        '  EINZELWERT-TEST, nachträglich eingefügt. Bitte die nachstehende Testzeile
        '  auskommentieren, falls nicht bestimmte Parameter getestet werden sollen!
        'F!=Ellpci( 0.9, 0.9 )  :clearclip:putclip str$(Ellpci_E!):waitinput
        '  Ergebnis in der Zwischenablage 17-stellig verfügbar, um Präzision zu prüfen.
        '*************************************************************************
        '*************************************************************************
        'Ausgabe:
        print " Phi = ";format$("00.0",&Loop);" °";
        print tab(15);"F(phi,k) = ";format$("%g",f!),
        print tab(45);"E(phi,k) = ";format$("%g",ellpci_e!)
        '*************************************************************************

    endwhile

Endwhile

waitinput
END
'{ VERIFIKATIONSTEIL
'*************************************************************************
' Präzisionsausgabe von E(phi,k) online verfügbar auf:
' https://keisan.casio.com/exec/system/1244989948
' Weitere Kontrolle via Maxima-Funktion elliptic_e(phi(=a),m(=k^2))
'*************************************************************************
'E(x,k) = Ellipci_E!
'-----------------------------------
'0.0, 0.0 = 0
'0.0, 0.1 = 0
'0.0, 0.4 = 0
'0.0, 0.9 = 0
'*** OK ! ***
'-----------------------------------
'Keisan 0.6,0.4   = 0.5945980735117437678766
'Maxima 0.6,0.4^2 = 0.59459807351174
'Ellpci(0.6,0.4)  = 0.59459807353087835
'*** OK ! ***
'-----------------------------------
'E(0.9,0.9) =
'Keisan(phi[rad],k):          0.808143018647826110606
'Maxima: elliptic_e(phi,k*k): 0.80814301864783
'Ellipci:                     0.80818524219652310
'*** MATHEMATISCH NICHT OK !!! (Technisch ausreichend) ***
'---------------------------------------------------------
'Ellipci(0.9,sqrt(0.9))           = 0.79659968020564376
'Maxima's_elliptic_e(0.9,0.9=k^2) = 0.79659968020564
'Keisan_E(0.9,Sqrt(0.9))          = 0.796599680205643585782
'*** OK ! ***
'Hinweis: Sqrt(0.9)=0.94868329805051379959966806332982
'-----------------------------------
'1.0, 0.0 = 1.570796326794896619231
'= Pi()/2 = 1.5707963267948966192313216916398
'*** OK ! ***
'}*********************************************************

Proc Ellpci :Parameters a!,k!

    ' a := Phi_max (sometimes called 'Amplitude')
    ' k := Modulus (In der Literatur wird auch m verwendet, mit m = k^2)
    ' F! = integral of first kind
    ' E! = integral of second kind
    ' Algorithm: Landen transformation with the DiDonato recurrence
    ' in specific ranges of the argument
    ' Proc uses external variable ellpci_e! for output of E!
    declare s! ,c!
    declare s0!,c0!, p0!,q0!
    declare s1!,c1!
    declare p1!,q1!, p2!,q2!
    declare s2!,s3!,s4!,s7!
    declare r0!,r1!,r2!, r7!
    declare n0!,n1!,n2!
    declare t0!,t1!,t2!,t3!,t7!
    declare d0!,d1!,d2!,d3!,d4!,d7!
    declare i1!,i2!,j1!,j2!,k1!,k2!, m7!,n7!
    declare L1!,L2!
    p2!=k!*k!
    q2!=(1-k!)*(1+k!)
    s!=sin(a!)
    c!=cos(a!)
    s0!=s!*s!
    c0!=c!*c!
    p0!=abs(k!*s!)
    q0!=abs(k!*c!)
    case p0!>=sqrt(.5):goto "ellpci_44"
    r1!=a!
    r2!=1
    s1!=0
    s2!=0
    n1!=1
    n2!=2
    t0!=a!*s0!
    t7!=s!*c!
    goto "ellpci_33"
    ellpci_29:
    n1!=n2!+1
    n2!=n1!+1
    t0!=s0!*t0!
    t7!=s0!*t7!
    ellpci_33:
    r0!=s1!
    r1!=(n1!*r1!-t7!)/n2!
    r2!=p2!*r2!/n2!
    s2!=s2!+r1!*r2!
    r2!=n1!*r2!
    s1!=s1!+r1!*r2!
    case abs(t0!)<abs(r0!):goto "ellpci_41"
    case abs(s1!)>abs(r0!):goto "ellpci_29"
    ellpci_41:
    f!=a!+s1!
    ellpci_e!=a!-s2!
    return f!
    print "Never reach 43-error!":waitinput
    ellpci_44:
    d7!=(1-p0!)*(1+p0!)
    d0!=sqrt(d7!)
    i2!=1
    j2!=1
    k2!=0
    m7!=0
    n7!=0
    s1!=0
    s2!=0
    s3!=0
    s4!=0
    t3!=q0!*d0!
    n0!=2
    goto "ellpci_63"
    ellpci_58:
    i2!=i1!
    j2!=j1!
    k2!=k1!
    n0!=n0!+2
    t3!=d7!*t3!
    ellpci_63:
    n1!=(n0!-1)/n0!
    n2!=(n0!+1)/(n0!+2)
    i1!=n1!*i2!
    j1!=n1!*n1!*q2!*j2!
    k1!=k2!+2/(n0!*(n0!-1))
    r0!=t3!/n0!
    m7!=n2!*n2!*q2!*(m7!-r0!*i1!)
    n7!=n1!*n2!*q2!*(n7!-r0!*i2!)'=Orig
    d1!=j1!
    d2!=n2!*j1!
    d3!=m7!-j1!*k1!
    d4!=n7!-n1!*q2!*j2!*k1! + q2!*j2!/(n0!*n0!)
    r0!=s3!
    s1!=s1!+d1!
    s2!=s2!+d2!
    s3!=s3!+d3!
    s4!=s4!+d4!
    case s3!<r0!:goto "ellpci_58"
    t0!=d0!+q0!
    r7!=t0!/4
    L1!= -1*ln(r7!)
    t7!=1+p0!
    s7!=t7!/2
    L2!= ln(s7!)
    'k=1:Kreis, daher Arc=Kreisbogen[rad] 'Meine eigene Änderung!

    if L2!=0:ellpci_e!=sin(a!):beep:return Ln(tan(a!/2+Pi()/4)):End:endif

        t1!=(1+s1!)*L1!+q0!/d0!*L2!
        t2!=(.5+s2!)*q2!*L1!+1-q0!/d0!*(1-p0!)
        f!=t1!+s3!
        ellpci_e!=t2!+s4!
        case a!>=0:goto "ellpci_94"
        f!= -1*f!
        ellpci_e!= -1*ellpci_e!
        ellpci_94:
        return f!

    endproc