VALIT, signaalverwerking voor Laser Doppler Anemometry
Citation for published version (APA):
Hommerson, G., & Koene, F. G. H. (1987). VALIT, signaalverwerking voor Laser Doppler Anemometry. Technische Universiteit Eindhoven.
Document status and date: Gepubliceerd: 01/01/1987 Document Version:
Uitgevers PDF, ook bekend als Version of Record Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at: openaccess@tue.nl
providing details and we will investigate your claim.
voor Laser Doppier Anemometry
G. Hommersom F.G.H. Koene
jan 87
Inhoudsopgave
Hoofdstuk 1 Inleiding
hoofdstuk 2 Matrixinversie
2.1 Computerkeuze voor de matrixinversie 2.2 Het inversieprogramma op de PRIME Hoofdstuk 3 Versturen van de geinverteerde matrix
van de PRIME naar de Apple
3.1 Versturen van de matrix van PRIME naar BURROUGHS
3.2 Versturen van de matrix van BURROUGHS naar Apple
3.3 Gebruiksklaar maken van de geinverteerde matrix
Hoofdstuk 4 Het rekenprogramma op de Apple 4. 1 Het hoofdprogramma (MAIN) 4.2 Procedure BOOT 4.3 Procedure REAOMATRI X 4.4 Procedure SETPARAMETER 4.5 Procedure INTAKE 4.6 Procedure MULTIPLYCORR 4.7 Procedure FIT 4.8 Procedures WRITETOFLOPPY en WRITETOPRINTER 4.9 Procedure ERRORHANDLER 4.10 Andere gebruikte procedures 4.11 Indices in het rekenprogramma
Hoofdstuk 5 5. 1 5.2 5.3 5.4 5.5
Verbinding tussen correlator en Apple Beschrijving van de gebruikte technieken Hardware aanpassing van de interface voor gebruik met de Apple
Correlator besturingssoftware Floating point berekeningen
Beschrijving van de gebruikte (externe) procedures
Programmalistings
Het matrixinversieprogramma MINV.F77 Het conversieprogramma FORMATCONVERSION Het verwerkingsprogramma VALIT
De correlatorbesturingsprogramma's
Het Assembler vermenigvuldigingsprogramma MATRIXMULTIPLY
In lit 1 is de signaalverwerking van autocorrelatiefuncties m.b.v. lineaire transformatie behandeld. Er wordt in deze bijlage ingegaan op de implementatie van die methode op de in de opstelling gebruikte Apple 2e.
Daarnaast wordt informatie gegeven over de aard van de communicatie tussen de Apple en de overige meetapparatuur en de realisatie daarvan.
Hoofdstuk 2. Matrixinversie
2.1 Computerkeuze voor het inverteren.
Het inverteren van de transformatiematrix (lit 1) is niet op de Apple PC gebeurt maar op een mainframe, de PRIME. Hiervoor waren verschillende redenen aanwezig.
1. Grotere geheugenruimte op de PRIME. De Apple is uitgerust met een geheugenuitbreidingskaart, waardoor het totale geheugen 128 kbyte is. Hiervan is een gedeelte nodig voor de systeembesturing. De
transformatiematrix en de geinverteerde daarvan benodigen aan opslag alléén al 71 kbyte (2*94*94*8 byte). Met inversieprogramma en
werkgeheugen was de totaal benodigde geheugenruimte wellicht te groot geworden.
2. Aanwezigheid van een standaardinversieprocedure in de NAG-bibliotheek van de PRIME. Deze bevat ook foutmeldingen voor slecht
geconditioneerde matrices.
3. Grotere precisie op de PRIME (64 bits tegenover 32 bits op de Apple). Dat is vooral van belang bij slecht geconditioneerde matrices. Bij het inverteren worden dan twee vrijwel gelijke getallen van elkaar afgetrokken.
4. Grotere snelheid. Op de PRIME duurt de inversie ongeveer 40 seconden aan rekentijd, de verwachte duur op de Apple enkele uren. Dat kan een probleem zijn bij het ontwikkelen en testen van de inversieprocedure. Nadeel van het gebruik van de PRIME is de noodzaak van het oversturen van de geinverteerde matrix naar de Apple, en het gebruiksklaar maken daarvan (hoofdstuk 3).
2.2 Het inversieprogramma op de PRIME.
Een flowschema van het in FORTRAN geschreven programma (MINV.F77) is geschetst in fig. 2.1. Voor verdere bijzonderheden wordt verwezen naar de van commentaar voorziene programmatext.
declaraties lezen invoer-parameters printen invoer-parameters nee correctie schaalparameters RJKG en TAUFKG 1 genereren matrix aanpassen indices matrix voor inversie matrixinversie (F01AAF) opslag geinverteerde matrix (stop)
Compileren en linken aan de bibliotheek gebeurt met CPL (Control Process Language) files, en wel door het commando: "R_COMLINKMINV.CPL". Uitvoeren van het programma met het commando: "R_RENMINV.CPL_INVER", waarbij voor de naam INVER ook een andere naam gekozen mag worden (bijv. INV20/4). De geinverteerde matrix wordt vervolgens geformatteerd
opgeslagen in de 4 files: INVER1.DATA t/m INVER4.DATA. Dit om het langdurig filetransport in stukken te hakken. Tenslotte wordt de geinverteerde matrix ook nog ongeformatteerd opgeslagen voor het niet langer gebruikte testprogramma SPELEN.F77.
Hoofdstuk 3 Versturen van de geinverteerde matrix van PRIME naar Apple.
3.1 Versturen van de matrix van PRIME naar BURROUGHS.
Omdat de PRIME niet rechtstreeks met de Apple kan communiceren, worden de files eerst verstuurd naar de BURROUGHS computer. Stel dat de te versturen files de namen INVER1.DATA, INVER2.DATA, INVER3.DATA en INVER4.DATA hebben (par. 2.2). Ze staan bijvoorbeeld op het kladgeheugen POOL>FRANS. De verbinding tussen PRIME en BURROUGHS wordt tot stand gebracht via een met het THE-net verbonden terminal. Eerst inloggen op de PRIME door de commando's:
- CALL_400 - LOGIN
- usercode (PRIME) - password (PRIME)
(verzoek aansluiting met de PRIME via het THE-net)
(vraag de PRIME om in te loggen)
Om het filetransport te beginnen de volgende commando's: - BURR
- usercode_password POOL (of b.v. USER5) - F - <return> (2x) - PB - POOL>FRANS>INVER1.DATA - INVER1 - <return> (4x)
- c
(aanroep van het filetransportprogramma) (van de BURROUGHS, om later in te loggen) (geheugenruimte op de BURROUGHS waar je
files terechtkomen) (kies voor Filetransport)
(als antwoord op vragen zoals: "bent U zeker van wat U intypt?")
(richting filetransport van PRIME naar BURROUGHS)
(naam van de eerste PRIME-file) (naam eerste BURROUGHS-file)
(default antwoord op die vragen is goed) (Copy; filetransport begint)
Inloggen op de BURROUGHS gebeurt nu automatisch waarna het
filetransport begint. Als de eerste file verstuurd is volgt dezelfde procedure, te beginnen met het commando •p•. Voor een 94x94 matrix hebben de vier PRIME files resp. 434, 452, 433 en 452 regels. Het versturen duurt ongeveer 4 minuten per file.
Als alle vier de files verstuurd zijn volgen een paar commando's ter afsluiting.
- Q - LO
(Quit, automatisch uitloggen op de BURROUGHS volgt)
(uitloggen op de PRIME)
Tenslotte een waarschuwing. Als op de BURROUGHS bij gebrek aan eigen geheugenruimte voor POOL gekozen is, moeten de files vóór zes uur 's avonds naar de Apple gestuurd zijn (en dat duurt ongeveer twee uur) omdat POOL dan gewist wordt.
3.2 Versturen van de matrix van BURROUGHS naar Apple
Voor het gestandaardiseerd versturen van files van BURROUGHS naar Apple bestaat het datatransportprogramma KERMIT. Helaas bleek dat niet te werken, zodat een minder modern transportprogramma FTP116 gebruikt is. Dit programma heeft als nadeel dat het werkt met een baud-rate
(snelheid van versturen van gegevens) van 300, wat niet razendsnel te noemen is. Voor het filetransport moet gebruik gemaakt worden van de Apple-A in het rekencentrum, die op het THE-net is aangesloten. De verbinding tussen die Apple en de BURROUGHS heeft gewoonlijk een baud-rate van 1200. Dat kan veranderd worden m.b.v. het programma KERMIT, dat wel met een baud-rate van 1200 werkt.
Start nu de Apple met in drive 1 de schijf KERMIT en in drive 2 een lege schijf. Als een 94x94 matrix wordt geinverteerd is voor elk van de 4 files 73 geheugenblokken nodig, wat dus niet past op één enkele schijf
(278 blokken). Na het versturen van drie files dus niet vergeten de schijf in drive 2 te verwisselen. Probeer de volgende commando's:
- x
- KERMIT - CONNECT - STA - BAUD_300 - <contr>-] C - EXIT- x
- FTP116 - STA (eXecutehet programma KERMIT)
(vraag verbinding met het THE-net) (geeft een overzicht van een stel
parameters, waarvan de baud-rate op 1200 moet staan; eventueel het commando
herhalen)
(verander de baud-rate. Als het commando niet begrepen wordt nog eens herhalen. Het is gelukt als een onzin antwoord komt in de trant van x•xx.)
(verbreek verbinding met THE-net) (stop het KERMIT programma)
(eXecute
het programma FTP116)
(baud-rate moet nu 300 zijn, schrik niet van dubbele tekens op het scherm)
- CALL_100 - usercode - password - FAMILY_DISK=POOL_ OTHERWISE_DISK - FILE - <esc> - A - #5:INVER1.TEXT - INVER1
(vraag verbinding met de BURROUGHS, even wachten tot je mag inloggen)
(jawel, van de BURROUGHS) (idem)
(naar de POOL geheugenruimte) (kijk of de files er staan) (begin het filetransport)
(transport in de richting van de Apple) (Apple filenaam; aanhangsel ".DATA" MAG
NIET)
(BURROUGHS filenaam)
Het versturen begint nu, regel voor regel.Na 20 regels (20 puntjes op het scherm) volgt een controle. Bij een fout (zelden) volgt een •E• op het scherm en worden de 20 regels opnieuw verstuurd, anders volgt een
•B•. Neem een goed boek mee, want het duurt ongeveer 30 minuten per file. Voor het versturen van de volgende file dezelfde procedure, te beginnen met <esc>.
Na afloop de troep opruimen: - BYE
- BAUD_1200 - <contr>-X
(uitloggen BURROUGHS)
(baud-rate terug naar 1200) (stop het programma FTP116)
De Apple is nog een beetje overstuur van al die commando's, waardoor eerst het <contr>-reset commando gegeven moet worden om het programma KERMIT weer te kunnen gebruiken (b.v. om te zien of de baud-rate weer terug op 1200 staat).
3.3 Gebruiksklaar maken van de geinverteerde matrix.
De 4 files zijn als karacters verstuurd, terwijl we ze graag
ongeformateerd op schijf willen hebben. Dat verbetert de leessnelheid door het rekenprogramma met een factor 30. Een en ander gebeurt met het programma CON, waarvan een flowschema in fig. 3.2 te vinden is.
Voor het gemak wordt in elk geheugenblok één matrixrij (94 real getallen) weggeschreven. Omdat een blok 128 real getallen kan bevatten, zijn de overige 34 plaatsen met nullen gevuld. Verder zijn de filenamen INVER1.DATA t/m INVER4.DATA default ingebouwd. Tenslotte bestaat het programma uit een lees- en een schrijfprocedure. Dit is gedaan om de foutmelding PROCEDURE TOO LONG te vermijden. Ook als er in het hele programma geen enkele procedure voorkomt (afgezien van standaard READ-en WRITE-procedures) kan die foutmelding optredREAD-en. Dat duidt er dan op dat teveel of te lange statements achter elkaar gebruikt worden.
files invoer invoer naam ongeformateerde opslag van één matrixrij
Fig. 3.2 Flowschema van het programma voor het ongeformateerd opslaan van de geinverteerde matrix
Hoofdstuk 4 Het rekenprogramma VALIT op de Apple. De gebruikte microcomputer is een Apple 2e met een
geheugenuitbreidingskaart, waardoor het totale geheugen 128 kbyte wordt. Uitgebreide documentatie is te vinden in (lit 2) en (lit 3). Voor de duidelijkheid volgen hier nog enkele bijzonderheden.
Van het rekenprogramma VALIT (Verwerking van Autocorrelaties m.b.v. Lineaire Transformatie) is een aantal representaties aanwezig.
Allereerst de textfile (achtervoegsel H.TEXT"), de gecompileerde versie daarvan (achtervoegsel ".CODE") en de aan de bibliotheken en externe
(Assembler) procedures gelinkte versie (ook achtervoegsel ".CODE"). Het programma bestaat in feite uit het achter elkaar aanroepen van procedures (par 4.1). Aangezien de benodigde ruimte voor de opslag van de textfile groter werd dan het maximale aantal van 50 blokken, zijn enkele procedures in een eigen bibliotheek gestopt, UNITPROC genoemd en voortaan simpelweg aangeduid met bibliotheek. Deze dient niet verward te worden met de reeds aanwezige standaardbibliotheek.
Voor het gebruik van de standaardbibliotheek moet aan het begin van het programma de statement •usES TRANSCEND" opgenomen zijn. Deze
statement mag niet 66k nog eens in de bibliotheek gebruikt worden. Procedures, waarin functies zoals EXP, COS of SQRT etc worden gebruikt kunnen daarom niet in de bibliotheek ondergebracht worden.
Alle TYPE en FILE declaraties moeten in de bibliotheek gebeuren. Alle variabelen zijn daarintegen (met uitleg) in het hoofdprogramma
gedeclareerd.
Een verschil tussen procedures in de bibliotheek en het
hoofdprogramma is dat van de eersten de procedureheading (met daarin de declaraties van in de procedure gebruikte globale variabelen) apart in het begin gedeclareerd moet zijn.
Tenslotte wordt een aantal externe (Assembler) procedures gebruikt. Dit zijn de procedures MATRIXMULTIPLY, voor het vermenigvuldigen van correlogram en matrix, en een aantal procedures voor de besturing van de MALVERN K7023 correlator: BOOTK7023, OPERATEK7023 en READK7023.
4.1 Het hoofdprogramma (MAIN)
Een flowschema van het hoofdprogramma wordt geschetst in fig. 4.1
BOOT READ-MATRIX commandoregel PARAM. SET nee ja IN-TAKE print nee PLY
Fig. 4.1. Flowschema van het op de Apple gebruikte verwerkingsprogramma.
Het programma is horizontaal opgebouwd, d.w.z. dat na het uitvoeren van elke opdracht (in de vorm van een procedure) wordt teruggesprongen naar de hoofdcommandoregel voor een nieuw commando. Als in een procedure een fout optreedt, wordt die zoveel mogelijk gespecificeerd.
Indien het uitvoeren van de opdracht foutloos gebeurd is krijgt de variable DEFAULTNEXT (een character) een nieuwe waarde. Die hoort bij de opdracht die bij normaal gebruik op de juist uitgevoerde volgt. Voor het uitvoeren van die opdracht hoeft slechts de <return>-toets ingedrukt te worden. Het door DEFAULTNEXT bepaalde flowschema is geschetst in fig. 4.2. Ten alle tijde kan echter zelf het volgende commando bepaald worden. READ-MATRIX PARAM. SET foutspec. ~---- foutspec.1----nee
Het printen van de hoofdcommandoregel en het verkrijgen van een nieuw commando gebeurt overigens in een aparte procedure GETVALIDCOMMAND.
4.2 Procedure BOOT
Deze procedure activeert de interface die de communicatie tussen de correlator en de Apple verzorgt. De const INIT (zie de declaratie in de bibliotheek) dient niet veranderd te worden.
4.3 Procedure READMATRIX
Een flowschema van deze procedure is te vinden in fig. 4.3. Naast de geinverteerde matrix worden enkele belangrijke parameters van de schijf gelezen, zoals het aantal kanalen van de correlator, het aantal fringes, de fringevisibility en of de matrix al dan niet geschaald is. Deze
parameters worden in het eerste blok gelezen. De parameters JMEANS(UARE) en TAUMEANS(QUARE) komen overeen met de variabelen RJKG en TAUFKG in het inversieprogramma. Bij het lezen van een blok m.b.v. de
standaardprocedure BLOCKREAD worden de gelezen matrixelementen (èèn rij) in het bufferarray REALBUFFER geplaatst. Voor het lezen van een nieuw blok moeten de waarden van die matrixrij op de goede plaats in het array INVMATRIX[I,J] gezet zijn
lees matrixnaam lees eerste blok nee lees volgend blok
Fig. 4.3. Flowschema van de procedure READMATRIX.
4.4 Procedure SETPARAMETER.
2
ja
Het flowschema van deze procedure is te vinden in fig. 4.4. Voor de meting van belang zijnde grootheden (o.m. de hoek ~ tussen de
laserbundels en de op de correlator ingestelde sampletijd T worden in
s deze procedure ingevoerd. Verder worden ook de uitvoerfiles (naar
printer en disk) geopend. Om de al eerder genoemde foutmelding PROCEDURE TOO LONG (par. 3.3) te vermijden is dit in een aparte procedure
SETOUTFILE gebeurt. Eventueel al aanwezige files worden overigens eerst gesloten. Op die manier kan bij elke aanroep van SETPARAMETER in een nieuwe file geschreven worden.
printen grootheden lezen van waarden voor grootheden berekening snelheden u. 1
4"
4.5 Procedure INTAKE.
Een flowschema is geschetst in fig. 4.5.
(
begin)
..
~
printen ;' ammandoregelI
lezenI
~ commando geldig nee commando? ja I I Jstart stop reset lees reset, start (stop
correlator correlator correlator correlator en lees correlator
,
•'
~Fig. 4.5. Flowschema van de procedure INTAKE.
Net als in het hoofdprogramma is er een DEFAULT variabele aanwezig (par. 4.1), zodat het aanslaan van de <return>-toets voldoende is voor het innemen van een correlogram en het verlaten van de procedure.
Het laatste element van het array AUTOCORRELATION wordt gelijk aan ·nul gemaakt. Dat komt door de extra vergelijking p
1
=
0 (lit 1).Aan de constanten RESET, START etc worden integerwaarden toegekend. Die worden gebruikt in de commando's met dezelfde naam. De constante DIMCORR, die overeenkomt met de integer LENGTHCORR in de
procedureaanroep, is in de bibliotheek gedeclareerd.
4.6. Procedure MULTIPLYCORR.
Een flowschema van de gebruikte procedure wordt in fig. 4.6 geschetst. correctie autocorrelatie autocorrelatie verminderen met gemiddelde vermenigvuld. autocorrelatie met matrix nee correctie kansverdeling piek kansver-deling zoeken kansverdeling normaliseren op 1 nee
Fig. 4.6. Flowschema van de procedure MULTIPLYCORR.
Inherent aan de methode van lineaire transformatie zijn de oscillaties in het hoge-snelheidsgedeelte van de berekende
snelheidsverdeling. Het is dus zaak om de piek in de eerste 60 à 70 kanalen van de snelheidsverdeling te houden. Dat kan met het veranderen van het rastertoerental. Daarmee verandert de snelheid waarmee de
fringes bewegen en daarmee de relatieve snelheid van de stroming t.o.v. de fringes (lit 1).
Om tijdwinst te boeken worden daarom ook alléén de eerste 72 waarden berekend. Dat aantal (de constante VELNUMCHAN; declaratie in de
bibliotheek) kan eventueel nog veranderd worden. Ook vanwege die
oscillaties kan het gebied waarin naar de piek gezocht wordt ingesteld worden m.b.v. de constanten LOWBOUND en UPBOUND (declaraties in de procedure).
Voor de vermenigvuldiging wordt de Assembler-procedure MATRIXMULTIPLY gebruikt. Om de snelheid daarvan op te voeren wordt het correlogram verminderd met het gemiddelde daarvan. Dat gemiddelde moet wel een heel getal zijn. De Assembler-procedure kan namelijk geen getallen tussen 0 en 1 verwerken.
4.7 De procedure FIT.
De parameters van de snelheidsverdeling worden geinterpoleerd in het snelheidsgebied rond de piek. Links van die piek worden POINTSLEFT punten meegenomen, rechts ervan POINTSRIGHT, met een maximum van MAXNUM aan elke kant. Alle punten waarvan de kans groter is dan 0.1 (de
maximale kans was op 1 genormaliseerd) worden meegenomen, en aan elke kant van de piek èèn extra. Een heel smalle piek heeft dus altijd drie punten. Als de ondergrens LOW-1 (-1 vanwege het extra punt) buiten de arraygrens dreigt te vallen wordt LOW-1 gelijk gemaakt aan die
arraygrens FIRSTINDEX.
Vervolgens worden de eerste 4 momenten berekend volgens formule 2.28 (zie verslag), en daaruit de gemiddelde snelheid U en standaarddeviatie
a volgens formule 2.29. De skewness (scheefheid) en de curtosis u
(spitsheid) van de snelheidsverdeling worden berekend volgens:
3 3
skewness
=
(m3- 3m1m2+ 2m1)/au enDe berekende waarden voor deze twee grootheden zijn echter onjuist. Dat in nagegaan met de analytische correlogrammen, gebaseerd op gaussische snelheidsverdelingen. De skewness daarvan moet 0 zijn, de curtosis 3.
Uitvoer naar het scherm gebeurt in een aparte procedure WRITETOSCREEN om de beruchte PROCEDURE TOO LONG foutmelding te verkomen.
4.8 Procedures WRITETOFLOPPY EN WRITETOPRINTER.
Indien besloten wordt tot opslag wordt een coordinaat van de meting gevraagd. Samen met een aantal van gelezen en berekende grootheden wordt het correlogram opgeslagen in èèn geheugenblok m.b.v. de
standaardprocedure BLOCKWRITE.
Het aantal gelukte experimenten wordt bijgehouden door de variabele NUMOFEXP. Na elke 50 experimenten begint de printer een nieuwe bladzijde en print daar o.a. de datum en een kop. Van elke meting wordt vervolgens de coordinaat, de gemiddelde snelheid, de standaarddeviatie, de skewness en de curtosis afgedrukt.
4.9 Procedure ERRORHANDLER.
In de meeste procedures wordt de foutvariabele INTFAIL (een integer) bijgehouden. Als na afloop van een procedure die variabele de waarde nul heeft wordt via DEFAULTNEXT (par. 4.1) naar de volgende procedure
gegaan. Als er wel een fout opgetreden is, heeft INTFAIL de bij die fout horende waarde gekregen. Daardoor wordt de procedure ERRORHANDtER
aangeroepen die de foutspecificatie op het scherm weergeeft.
De gehele procedure is weer in drieen gedeeld vanwege de foutmelding (jawel) PROCEDURE TOO LONG.
4.10 Andere gebruikte procedures.
Er is geconstateerd dat het tekenprogramma TURTLEGRAPHICS niet werkte indien er in hetzelfde programma een leesopdracht van disk aanwezig is. Om figuren zoals de berekende snelheidsverdeling op zijn minst te kunnen schetsen is gebruik gemaakt van de mogelijke cursorposities op het
scherm (80 in horizontale en 24 in verticale richting). De
snelheidsverdeling wordt zo m.b.v. de procedure SKETCHVEL na het vermenigvuldigen van matrix en correlogram geschetst. De mogelijkheid bestaat ook het ingenomen correlogram voor die vermenigvuldiging op dezelfde manier te schetsen m.b.v. de procedure SKETCHAUTO.
Een drietal procedures is ondergebracht in het AUXILIAIRY commando. Allereerst is er de procedure GENERATECORR, die volgens formule 2.16
(verslag) een autocorrelatie gegenereert. Er worden twee e-machten berekend. Om underflow te vookomen wordt nagegaan of het argument van die e-machten groter is dan -50. Het invoeren van de parameters gebeurt in een aparte procedure PARAM.
Verder is er de procedure READCORR, die een eerder opgeslagen
correlogram van schijf leest. Een flowschema is geschetst in fig. 4.7. Na de laatste test wordt onafhankelijk van de uitkomst daarvan -naar het einde van de procedure gesprongen. In geval er een fout
opgetreden is krijgt INTFAIL een van nul afwijkende waarde en zal door de aanroep van de procedure ERRORHANDtER de fout gespecificeerd worden.
Tenslotte is er nog de procedure CONTRPAR, die desgewenst alle van schijf gelezen en in PARAMETERSET ingevoerde parameters laat zien.
lees filenaam nee ja~---~ lees nummer geheugenblok nee
4.11 Indices in het rekenprogramma.
In de onderstaande tabel wordt een overzicht gegeven van het bereik van een aantal indices
omschrijving bereik en overeenkomstige waarden indices
interne correlatorkanalen 0 1 2 3 4 5 6 7 8 99
index array CORRELATOR 1 2 3 4 5 6 7 8 9 100
aantal t vertaging 1 2 3 4 5 96
s
index array AUTOCORRELATION 1 2 93 94
1 2 93 94
index array INVMATRIX[i,j]
72 93 94
..
index array VELOCPROB 1 72
Tabel 4.1 Bereik van indices in het verwerkingsprogramma.
4.12 Veranderingen in het rekenprogramma.
Hier moet gewezen worden op een verschijnsel dat niet verklaard is kunnen worden. Na veranderingen moet het rekenprogramma opnieuw
gecompileerd en gelinkt worden. Zelfs als dat zonder foutmeldingen geschied is en er een analytisch correlogram (bijv. U
=
30, a=
5 m/s)u ingevoerd wordt, kan het voorkomen dat het programma een foutieve snelheidsverdeling berekent (fig. 4.8).
+
••
... +•
+ ~b
t +•
+d
•
a
+•
c
+ + f•
..
+•
••
t +•
•
....
+ +..
1' t..
+ t ++ t..
f' + + + t + ... + + t i' +••
+ t.
....
•
.....
+..
•
• + • t..
... ++++ +..
...
++++ t tt+ + ++Fig. 4.8 Berekende snelheidsverdelingen, gewenst (a) en foutief (b,c en d)
Mogelijk is de oorzaak hiervan dat het programma te groot geworden is en er ergens een verkeerd geheugen geadresseerd wordt. Opnieuw compileren
(inclusief bibliotheek) en linken kan wel een goed resultaat opleveren. Op de schijf BACKUP is steeds een versie aanwezig gehouden die geen van de bovenbeschreven problemen opleverde. GOOI DIE NIET WEG.
Het verdient aanbeveling om na een verandering van het programma ook de bibliotheek opnieuw te compileren. Dat moet vóór het eigenlijke programma gebeuren. Start de Apple met in drive 1 (te adresseren met "#4:" de schijf Apple1 en in drive 2 (adres "#5:") de schijf Apple2. Geef daarna de volgende commando's:
- c
- #4:UNITPROC - #5:ANALCORR.LIB.
- c
(Compile)
(staat op schijf Apple1 dus vertel de com-puter dat hij in drive 1 moet zoeken) (zo heet de gecompileerde bibliotheek;
vergeet de laatste punt niet) (Compile)
- #S:ANALCORR - $ - L - <return> - #S:ANALCORR.LIB. - M1 - K7023 - <return> (2x) - #S:VALIT
(staat op schijf Apple2 dus "#5:" ervoor) (de gecompileerde versie wordt dan onder
de naam ANALCORR.CODE op schijf Apple2 gezet)
(Link)
(linker opent automatisch ANALCORR.CODE) (bibliotheek eraan linken)
(externe procedures eraan; M2 mag ook maar die blijkt niet altijd te werken)
(nog meer externe procedures) (er hoeft verder niets meer aan) (de naam van de gelinkte versie; mag
natuurlijk ook een andere zijn) Het programma wordt gestart met de commando's:
- x
(eXecute)- #S:VALIT
Als het programma onder de naam "SYSTEM.STARTUP" op schijf aanwezig is, wordt het programma automatisch gestart na het aanzetten van de Apple.
Hoofdstuk 5 Verbinding tussen correlator en Apple
Een Malvern K7023V autocorrelator kan naast besturing vanaf het frontpaneel ook extern bediend worden. De meegeleverde
correlatorprocessor bezat deze mogelijkheid.
Het VALIT autocorrelatie verwerkingsprogramma dient toegang te hebben tot de inhoud van de geheugenkanalen van de autocorrelator. Daar de meegeleverde MOP 7023 computer en de Apple 2e geconstrueerd zijn rond dezelfde (6502) processor, is de originele interfacekaart gebruikt voor de communicatie tussen Apple en Malvern K7023. Het uitvoerprotocol van de correlator wordt beschreven in lit 7.
5.1 Beschrijving van de gebruikte technieken.
De controle in- en uitvoerchip (INTEL 8255A) van de interfacekaart wordt bestuurd vanuit de Apple door middel van de •direct memory access• techniek. De in- en uitgang van de 8255A-chip staan via
uitbreidingspoort nr 7 in directe verbinding met de adresbus van de Apple 2e. De 8255A-chip is geactiveerd indien èèn van de
uitbreidingspoortadressen (hexadecimale notatie $COFC t/m $COFF) worden geadresseerd door de Apple. Hierbij hebben de poortadressen de volgende functie:
$COFC: programmerstand INTEL 8255A, de gewenste configuratie van de in-en uitvoerchip kan wordin-en gekozin-en.
Bij de gekozen configuratue (mode 0, controle woord #3) wordt de functie van de andere geheugenlocaties:
$COFD: monitorpoort. waarbij de bits 0-3 gegevens vanuit de correlator doorgeven.
bit 0: correlator klaar: ja - 0 nee - 1 bit 1: data stabiel: ja - 1
bit 2: externe triggering (niet gebruikt) bit 3: correlator present (niet gebruikt)
bits 4-7 geven informatie van de Apple naar de correlator. De primaire functie is het coderen van de
uitvoerconfiguratie van de correlator, waarbij bit 4 als extra functie het commando "uitvoer gewenst" doorgeeft indien het van laag naar hoog geschakeld wordt.
$COFE: gegevenspoort, d.w.z. dat alle bits worden gebruikt om de inhoud van de correlator geheugenkanalen cijfer voor cijfer in ASCII-code door te geven aan de Apple.
$COFF: commandopoort: transmissie van bedieningscommando's van de Apple naar de correlator.
bit 0: laag - hoog - laag: correlator reset bit 1 : laag - hoog - laag: correlator start bit 2: laag - hoog - laag: correlator stop bit 3-5: niet in gebruik
bit 6: hoog: uitlezing gewenst
bit 7: laag: processor klaar voor uitlezing laag - hoog - laag: volgend byte
5.2 Hardware aanpassing van de interface voor gebruik met de Apple. De architectuur van de Apple 2e leende zich voor het voor het gebruik van de OMA techniek, die door Malvern was toegepast bij het ontwerp van de interface. Herprogrammering van de gebruikte geheugenlocaties voor de OMA was noodzakelijk. Hiervoor is de verbindingskaart in
uitbreidingspoort 17 van de Apple zodanig opgebouwd dat de interfacepoort de adressen $COFC t/m $COFF interpreteert als de oorspronkelijk in de MOP 7023 gereserveerde $AFFC t/m $AFFF.
5.3 Correlator besturingssoftware.
Het UCSD-Pascal ontbeert de mogelijkheid om geheugenadressen direct aan te spreken. Derhalve was het noodzakelijk om alle correlator
besturingshandelingen via machinetaal programma's te laten verlopen. De listings van de procedures geven in principe voldoende informatie over de functie en werkwijze van de procedures. Informatie over de
machinetaal kan gevonden worden in lit 5, over de systeemarchitectuur, het gebruik van uitbreidingspoorten in lit 6. De communicatie tussen Assembler en Pascal programma's wordt beschreven in lit 2 en lit 3.
5.4 Floating point berekeningen.
Een aanvulling op de in de vorige paragraaf genoemde literatuur is gewenst m.b.t. het format waarin "real" getallen worden opgeslagen: het
"floating point format•. Een real wordt opgeslagen in vier bytes, als geschetst in fig. 5.1. byte 0
[
I I I I I mantisse byte 3 l i l l i l I • exponent -il>f
teken Fig. 5.1. Representatie van een •real" getal in vier bytes.De exponent is in •two's complement" notatie geschreven (lit 5). Alle floating point berekeningen in de hier beschreven machinetaal procedures hanteren een iets ander format voor interne berekeningen. De mantisse staat in de eerste 3 bytes van het getal en de exponent in de vierde byte. Het teken wordt (indien noodzakelijk) elders bijgehouden. De stroomschema's van de gebruikte rekenhandelingsalgoritmen zijn weergegeven in fig. 5.2 t/m 5.4.
I
tel mantisse opl
l
1
mantissebits 1 positie lager Jl
zet carry in hoogste mantissebit
I
t
verhoog exponent met 1J
Fig 5.2a. Floating point optelling indien de exponenten gelijk zijn. verhoog laagste exponent met 1
mantissebits 1 plaats lager
1 in hoogste mantissebit
0 in hoogste mantissebit
verschuif mantisse-bits 1 plaats lager
vergelijk bytes mantisse
trek grootste mantisse van de kleinste af schuif mantisse 1 bit omhoog verlaag exponent met 1 ee antwoord 0
Fig. 5.3a. Floating point aftrekking als de exponenten gelijk zijn.
1 in hoogste mantissebit
verhoog laagste exponent met 1 mantisse 1 bit lager
0 in hoogste mantissebit
2 trek mantisse kleinste getal af nee schuif mantisse 1 bit omhoog verlaag exponent met 1
2
nee
converteer expon. f1 naar 'absoluut' format
schrijf f2 in antwoord teller:= lengte mantisse expon. antwoord:= expon. f1 + expon. f2 ja mantisse f2 1 bit lager
hoogste bit:=1 0 ingevoegd
2 mantisse f1 positie hoger ja nee tel mantisse f2 op bij mantisse antwoord
nee ja mantisse antwoord 1 bit lager expon. antwoord 1 hoger mantisse f2 1 bit lager
5.5 Beschrijving van de besturingsprocedures. BOOTK7023 (command, errorcode).
Als command de waarde 0 krijgt wordt de in- en uitvoerchip op de interfacekaart geprogrammeerd. Vervolgens wordt de juiste
uitvoerconfiguratiecode op de monitorpoort gezet, de commandopoort op 0 gezet, de correlator op "stop" gezet en maximaal 3 ~s gewacht op het "correlator klaar signaalu. Verloopt alles korrect, dan blijft de foutcode gelijk aan 0. Komt het •correlator klaaru signaal niet binnen de gestelde tijd dan heeft de foutcode de waarde 9. Indien command niet de waarde 0 heeft wordt de programmering van de I/0-chip achterwege gelaten.
OPERATEK7023 (command).
De waarde command dient zodanig gekozen te worden dat het bitpatroon direct op de commandepoort gezet kan worden om het bit •aan• te zetten. De procedure zet het bit zelf weer •uit". Eénduidige waarden voor
command zijn derhalve: 1 =reset (bit 1) 2 = start (bit 2) 4 =stop (bit 3).
READK7023 (correlogram[dim2], dim2, errorcode). De procedure READK7023
- bestuurt de correlator tijdens tijdens de uitleesprocedure, - leest de uitgang van de correlator,
- reconstrueert de geheugeninhoud en geeft deze in floating point format rechtstreeks door aan VALIT en
- genereert een foutmelding indien
- het "correlator klaar" signaal niet volgt na het "stop• signaal, - de wachttijd per kanaaluitlezing (0.26 s) wordt overschreden
(de foutcode wordt het kanaalnummer + 10) en
- de uitvoer niet past in de correlogramvector, en wel foutcode 8 als de vector te klein is en foutcode 7 als de vector te groot is. READK7023 wacht na aanroep tot het "correlator klaar" signaal wordt aangetroffen op de monitorpoort. De gebruiker kan de wachttijd verkorten
door een toets aan te slaan. De correlator krijgt dan een stopcommando vanuit KEYHIT. Volgt geen •correlator klaar• signaal dan geeft een tweede toetsaanslag de betreffende foutmelding op. Volgt wel een
•correlator klaar• signaal dan worden de vereiste commando's voor uitlezing op de monitor en controlepoorten gezet.
Het uitvoerprotocol van de correlator wordt beschreven in lit 7. De correlatoruitvoer is in ASCII-code waarbij de cijfers van het getal van hoog naar laag worden doorgegeven. Transformatie naar floating point format geschiedt volgens het stroomschema van fig. 5.5. Het
getransformeerde cijfer wordt opgeteld bij het resultaat, waarna het totaal met 10 wordt vermenigvuldigd (TIMETEN algoritme: a*10
=
(a*2*2 +a)*2 ) en het volgende getal wordt gelezen. Nadat alle cijfers zijn gelezen volgt de conversie naar UCSD-Pascal floating point format met
•t• teken en wordt het weggeschreven op het geheugenadres van het betreffende correlogramelement. Het programma stopt indien:
- het •end of file" teken •*• van de correlator wordt gelezen en - het aantal correlogramelementen overeenkomt met de voorgeschreven
dimensie (dim2).
Ontbreekt één van beide criteria volgt een foutmelding, met foutcode 7 indien het "end of file" teken gelezen is vóór dim2 bereikt is en met foutcode 8 in het tegenovergestelde geval (dim2 te klein).
expon. antwoord:=4 (131 in two's compl.)
mantisse antwoord:=O lees ASCII cijfer
schuif cijfer 4 bits omhoog schuif cijfer 1 bit omhoog exp. buffer:= exp. buffer -1 ja ja antwoord 0 wordt niet gebruikt rest cijfer = mantisse
MATRIXMULTIPLY (matrix, correlogram, snelheidsverdeling, dim1, dim2). Dit is een procedure die een matrix met een vector vermenigvuldigt. Zij is ontworpen voor vermenigvuldiging van de snelheidsverdeling, waarbij gebruik is gemaakt van het feit dat in het correlogram slechts gehele getallen voorkomen in floating point format. Dat houdt in dat er een verband bestaat tussen het aantal bits van de mantisse die informatie bevatten en de grootte van de exponent. Bij de vermenigvuldiging wordt dan ook alleen dát gedeelte van de mantisse bij de berekeningen
meegenomen. Hoe kleiner het getal, hoe minder informatiedragende mantissebits er zijn en dus hoe sneller de verwerking is. Dat is de reden waarom in de procedure MULTIPLYCORR het gehele gedeelte van het correlogramgemiddelde wordt bepaald en afgetrokken van het gemeten correlogram. Door deze handelwijze wordt de verwerkingstijd met ca 20\ verkort.
Verwijzigingen lit 1 lit2 lit3 lit4 lit 5 lit 6 lit7
Ontwikkeling van een rekenprogramma voor de analyse van autocorrelatiefuncties m.b.v. lineaire transformatie F. Koene
afstudeerverslag TUE, januari 1987, WV 157-029 Apple Pascal
Language Reference manual (met 128k uitbreiding) Apple Pascal
Operating System Reference Manual (met 128k uitbreiding) Studie van de swirlbeweging bij directinjectie dieselmotoren d.m.v. Laser Doppier Anemometrie
P. Snauwaert
promotieverslag RU Gent, 1984
Programming the Apple in Assembly language Rodnay Zaks
Berkeley, 1985
Apple reference manuel
Operating and Installation manuel, type K7023 Malvern digital correlator
INTEGER I,J,NK,NM,IA,IAINV,IFAIL,IIN,IOUT,IBEGIN,SCHAAL, ALFA
C I, J TELLERS
C NK AANTAL KANALEN CORRELATOR C NM DIMENSIE MATRIX
C IA,IAINV NODIG VOOR INVERSIEROUTINE (ZIE HANDLEIDING F01AAF) C IFAIL INTEGER DIE OPTREDENDE FOUTEN SPECIFICEERT
C IIN INVOER VAN SCHERM VOOR IIN=1 C JOUT UITVOER NAAR SCHERM VOOR IOUT=1
C !BEGIN KANAALNUMMER VAN HET EERSTE MEEGENOMEN CORRELOGRAMELEMENT C SCHAAL GESCHAALDE MATRIX: SCHAAL=1, ANDERS 0
C ALFA FRINGEVISIBILITY (OOK MODULATIEDIEPTEl REAL•B Pl,NF,TAUFKG,RNK,RI,RJ,RJKG,A,AINV,WKSPCE,
* TAUF
C PI 3. 14 ..
C NF AANTAL FRINGES IN HET HALVE MEETVOLUME
C TAUFKG TAU•F KWADRAAT GEMIDD. CORRECTIE VOOR GESCHAALDE MATRIX C RNK REAL GETAL VOOR HET AANTAL KANALEN
C RI,RJ REAL GETALLEN VOOR DE TELLERS
C RJKG REAL GETAL VOOR J KWADRAAT GEMIDD. ZIE TAUFKG
C A MATRIX
C AINV GEINVERTEERDE MATRIX
C WKSPCE WORKSPACE VOOR DE INVERSIE ROUTINE C TAUF TAU * FREQUENTIE
REAL RSCHAAL,RIBEGIN,RALFA COMMON/EEN/A(1:129,1:129) COMMON/TWEE/AINV(1:129,1:129) COMMON/DRIE/WKSPCE(1:129),TAUF(1:129) $INSERT SYSCOM>ASKEYS IIN=1 IOUT=1 WRITE(IOUT,5)
5 FORMAT(1X, 'VOER IN ALS INTEGERS: '/
1X, 'AANTAL KANALEN CORRELATOR, MAXIMAAL 128') READ(IIN,*)NK WRITE(IOUT, 10) 10 FORMAT(1X,'BEGINKANAALNR. CORRELATOR') READ(IIN,*)IBEGIN WRITE(IOUT,20) 20 FORMAT(1X,'MODULATIEDIEPTE (\)') READ(IIN, *)ALFA WRITE(IOUT,30)
30 FORMAT(1X, 'VOER IN ALS DOUBLE:'/
1X, 'AANTAL FRINGES IN BUNDELSTRAAL (BIJV. 13.5DO)') READ(IIN, * )NF
WRITE(IOUT, 35)
35 FORMAT ( 1X, 'ALS JE WILT SCHALEN, TYPE DAN 1' ) READ(IIN,*)SCHAAL
WRITE(IOUT,40)NK,NF,ALFA,IBEGIN
40 FORMAT(1X,'AANTAL KANALEN CORR.:',3X,I3/ 1X,'AANTAL FRINGES:',BX,F4.1/ c c c c c c c c 50 60 90 100 150 160 170 180 c c NM=NK+1-(IBEGIN-1) RNK=DBLE(FLOAT(NK)) • • IF (SCHAAL.EQ.1) THEN WRITE(IOUT,50)
FORMAT(1X, 'MATRIX GESCHAALD') RJKG=(1.0DO+RNK*RNK)/2.0DO
TAUFKG=(1.0DO/(RNK•RNK)+0.25D0)/2.0DO ELSE
WRITE (I OUT, 60)
FORMAT(1X, 'MATRIX NIET GESCHAALD') RJKG=O.ODO
TAUFKG=O.ODO ENDIF
GENEREREN VAN DE MATRIX A DO 100 I=IBEGIN,NK RI=DBLE(FLOAT(I)) TAUF(I)=((RI+1.0DO)*RNK-2.0DO*RI)/(2.0DO•RNK*(RNK-1.0DO)) DO 90 J=IBEGIN,NK RJ=DBLE(FLOAT(J)) A(J,I)=(1.0DO+(ALFA/100.)*DCOS(2.0DO*PI*TAUF(I)*RJ))* DEXP(-(TAUF(I)*TAUF(I)-TAUFKG)*(RJ*RJ-RJKG)/ (NF•NF)) CONTINUE CONTINUE DO 150 J=IBEGIN,NK RJ=DBLE(FLOAT(J)) A(J,NK+1)=DEXP(TAUFKG*(RJ*RJ-RJKG)/(NF*NF)) CONTINUE A(NK+1,IBEGIN)=DEXP((TAUFKG*((RNK+1.0DO)*(RNK+1.0DO)-RJKG)+ RJKG/(RNK*RNK))/(NF*NF)) DO 160 I=IBEGIN+1,NK+1 A(NK+1,I)=O.ODO CONTINUE
AANPASSEN VAN DE INDICES OMDAT DE INVERTEER-PROCEDURES VERONDERSTELLEN DAT DE EERSTE INDEX VAN DE MATRIX (A) 1 IS.
DO 180 I=1 ,NM D0170J=1,NM A(J,Il=A(J+IBEGIN-1,I+IBEGIN-1) CONTINUE CONTINUE IA=129 IAINV=129 IFAIL=O CALL F01AAF(A,IA,NM,AINV,IAINV,WKSPCE,IFAIL)
220 FORKAT(1X,'INVERTEREN IS GELUKT,'/
1X,'GEINVERTEERDE MATRIX IN OPGEGEVEN FILE(S)')
c
C ONGEFORMATTEERD WEGSCHRIJVEN OM TE GEBRUIKEN IN DE PRIME C PROGRAMMA'S c c WRITE(5)SCHAAL,IBEGIN,NK,ALFA, SNGL(TAUFKG),SNGL(RJKG),SNGL(NF), ((SNGL(AINV(I,J)),I=1,NM),J=1,NM) C GEFORMATTEERD WEGSCHRIJVEN (IN STUKKEN) VOOR DATA-C TRANSPORT MET FTP116 c RSCHAAL=FLOAT(SCHAAL) RIBEGIN=FLOAT(IBEGIN) RALFA=FLOAT(ALFA) WRITE(6,240)RSCHAAL,RIBEGIN,SNGL(RNK),RALFA, SNGL(TAUFKG),SNGL(RJKG),SNGL(NF), ((SNGL(AINV(I,J)),I=1,NM),J=1,NM/4) WRITE(7,240)((SNGL(AINV(I,J)),I=1,NM),J=NM/4+1,NM/2) WRITE(8,240)((SNGL(AINV(I,J)),I=1,NM),J=NM/2+1,3*NM/4) WRITE(9,240)((SNGL(AINV(I,J)),I=1,NM),J=3*NM/4+1,NM) 240 FORMAT((5(E13.6,1X))) END IF CALL EXIT END
(*$G+*) const firstindex lastindex type matrix var invmtrx jms numfringes ordermatrix rfrstchan riscale rmdepth rnumchan taums
lallows goto statement)
= 1; =97;
=array[firstindex .. lastindex,
firstindex .. lastindex] of real; :matrix; :real; :real; :integer; :real; :real; :real; : real; :real; {the parameters to be lread from the input !files and written to {the output file
procedure readfiles(var invmatrix :matrix;
var dummystring endoffile fileint firstchan i intfail infile infilename j numofchan numoffiles ready space var var var var var var var var jmeansq numoffringes orderotmatrix rfirstchan rintscale rmoddepth rnumofchan taumeansq :string; :char; :integer; :integer; :integer; :integer; :interactive; :string; :integer; :integer; :integer; :char; :char; :real; :real; :integer; :real; :real; :real; :real; :real);
IEOF becomes true AFTER reading the endoffile characterl lloop index}
{integer value of rfirstchanl lloop index}
{specifies type of error) !input file(s)l
{loop index)
{number of input files)
lto read the last space in the file begin lof procedure readfiles)
writeln('enter the number of files to be read'); readln(numoffiles);
(*$!-*) {turn off automatic l/0 checking)
end;
dummystring:=concat('IS:INVER',dummystring,'.TEXT'); repeat
writeln('enter the name of the next file to be read, '); writeln('type <sp> for ',dummystring);
readln(infilename); if infilename=' ' then
infilename:=dummystring;
writeln('put in drive 2 the floppy containing ',infilename); writeln('wheri ready, press return');
readln(ready); reset(infile,infilename); intfail:=ioresult; if intfail<>O then begin end if intfail=10 then
writeln(infilename,' not on specified volume') el se
writeln('I/0 error ',intfail) until intfail=O; writeln(' ... 1 ... 2 ... 3'); if fileint=1 then read(infile,rintscale,rfirstchan,rnumofchan,rmoddepth, taumeansq,jmeansq,numoffringes); firstchan:=trunc(rfirstchan); numofchan:=trunc(rnumofchan); orderofmatrix:=numofchan+1-(firstchan-1);
for j:=((fileint-1)*orderofmatrix) div numoffiles +1 to (fileint*orderofmatrix) div numoffiles do begin
end;
wri te ( ' . ' ) ;
for i:=1 to orderotmatrix do read(infile,invmatrix[i,j)) writeln;
read(infile,space,endoffile); if not eof(infile) then
writeln('end of file not reached'); intfail:=ioresult;
if intfail<>O then
writeln('l/0 error ',intfail); close(infile);
note(40,20)
(*$!+*) {turn automatic errorchecking back on) end; lof procedure readfilesl
procedure writefile(var invmatrix var jmeansq var numoffringes var orderotmatrix :matrix; :real; :real; :integer;
var var rnumofchan taumeansq :real; :real); label 99; var blocknumber blockstransferred i intfail j outfile outfilename ready realbuffer :integer; :integer; :integer; :integer; :integer; :file; :string; :char; :array[1: 128] of real;
lnumber of transferred blockl lnumber of blocks transferredl {loop index}
{specifies type of error} {loop index}
{untyped output filel
{buffer array used in blockwritel begin lof procedure writefilel
I *SI-• l
repeat
writeln('enter name of output file'); readln{outfilename);
writeln('put in drive 2 a floppy with at least ',orderofmatrix+1, ' free blocks');
writeln('when ready press return'); readln(ready);
rewrite(outfile,outfilename); intfail:=ioresult;
if intfail<>O then ·
writeln('I/0 error, intfail= ',intfail) until intfail=O; realbuffer[1]:=rintscale; realbuffer[2):=rfirstchan; realbuffer[3):=rnumofchan; realbuffer[4]:=rmoddepth; realbuffer[5]:=taumeansq; realbuffer[6]:=jmeansq; realbuffer[7]:=numoffringes; for i:=8 to 128 do realbuffer[i]:=O.O; blockstranferred:=blockwrite(outfile,realbuffer,1,0); intfail:=ioresult;
if (intfail<>O) or (blockstransferred<1) then begin
end;
writeln('error in writing block 0'); goto 99;
for i:=1 to orderotmatrix do begin blocknumber:=i; for j:=1 to orderafmatrix do end; blockstransferred:=blockwrite(outfile,realbuffer,1,blocknumber); intfail:=ioresult;
if (intfail<>Ol or (blockstransferred<1) then begin
end
writeln('error in writing block ',blocknumber); goto 99
99: close{outfile,lock);
( *$I+•)
if intfail<>O then
writeln('I/0 error ',intfail) end; lof procedure writefilel
begin {of main program}
readfiles(invmtrx,jms,numfringes,ordermatrix,rfrstchan,riscale,rmdepth, rnumchan,taums);
writefile(invmtrx,jms,numfringes,ordermatrix,rfrstchan,riscale,rmdepth, rnumchan,taums);
(*$G+*) label 99; var airtem autocorr blcknumber boostpres callreadmatrix
lallows the use of the goto statement!
: real; I air temperature I
:corvector;lvector of autocorrelationl
:integer; lnumber of block read in readmatrixl
:real; lboostpressurel
:boolean; ltrue if procedure readmatrix has
been calledl
callsetparameter :boolean; ltrue if procedure setparameter has
been calledl continue curtos day defaultnext deltapres dimeasure fileid fitchoise frstchan freqgrating frinspacing ifail ipeak invmtrx jms lload meanvel mdepth month next noise numchan numexp numfringes numpoints order•atrix outname perpvel possvel rev rradius skewnes standev :boolean; :real; :integer; :char; :real; :char; :string; :char; :integer; :real; :real; :integer; :integer; :matrix; :real; :real; :real; :integer; :integer; :char; :integer; :integer; :integer; :real; :integer; leurtos is I lday of monthl
(default next commandl lpressure difference orificel
ldirection of measurement (rad, tang or axl !name of file containing
inverted matrix!
I'M':fit with moments; 'G' :gaussian curve fitl (number of first meaningful
autocorrelation channell lfrequency of gratingl lfringespacingl
(specifies type of error)
lpeakindex of velocity distributionl linverted 11atrixl
lcorrection factor for scaled matrix! lload of enginel
(mean velocityl
(percentage of modulation depthl lnext commandl
lnumber of channels of correlatorl lnumber of experimentsl
lnumber of fringes in beam radius) lnumber of points used to calculate
the fit of the velocity distributionl :integer; lorder of matrix!
:string; (name of disk output filet
:real; (velocity perpendicular to ut
:velvector; !vector of possible veloeities u, :real;
:real; :real; :real;
depending upon the value of taul lrpm of enginel
lcoordinatel
lstandarddeveation of velocity distributionl
trigcount ttau taums velprob" visib :integer; :real; :real; :vel vector; :integer;
ltrigger counter for operation of correlatorl (time lag of autocorrelationl
lcorrection factor for scaled matrix! (vector containing velocity probability distributionl
procedure MATRIXMULTIPLY(mat:matrix; var int:corvector; var out:velvector;
· dim1, dim2: integer); external;
procedure BOOTK7023(mode:integer; var ifail:integer);external; procedure OPERATEK7023(mode:integer);external;
procedure READK7023(var dataout:rawdata; length:integer; var ifail:integer);external;
var boostpressure :real;
var dayofmonth :integer;
var deltapressure :real;
var directionmeasure :char;
var firstchan :integer;
firstindex :integer;
var freqofgrating :real;
var fringespacing :real;
var intfail :integer;
var load :real;
var monthofyear :integer;
var numofchan :integer;
var numofexp :integer;
var passvelocity :velvector;
var rpm :real;
var tau :real;
var triggercount :integer;
velnumofchan :integer);
const
numoflines =16384; {number of radial lines on the gratingl
var angle const1 const2 i pi preshiftvelocity
:real; langle between laser beams in gradl
:real; :real; :integer;
:real; l3o1417ooool
:real; begin lof procedure setparameterl
writeln('enter the following parameters'); writeln;
writeln( 'day of month: o o o o o o o o o o o o o o o o o ');
writeln('month (as integer): o o o o o o o. 0 0 o o 0 o');
writeln('tau (in nanoseconds):o 0 o o 0 0 0 0 0 0 0 0 o');
writeln('frequency of grating (in Hz): o o 0 0 o o 0 o');
writeln('angle between laser beams (in grad): o');
lwriteln('trigger count:o o o o o o o o o o 0 0 o o 0 o');
writeln('rpm of engine (rpm): o o o o o o o o 0 o o 0 o');
writeln( 'boostpressure (mm Hg): o o o o o 0 o o o 0');
writeln('pressure diffo of orifice (mm Hg): 0 0 o 0 0 0 ' ) ;
writeln('temperature air (C):o o o o o o o 0 0 o o 0 o');
writeln('load (kilogram): o o o o o o o o o o o o o o o ');I
writeln('T(angential, R(adial or A(xial measuremento 0 ' ) ;
gotoxy(47,4); readln(dayofmonth); gotoxy(47,5); readln(monthofyear); gotoxy(47,6); readln(tau); if tau< 50 then intfail:=240; tau:=tau*1E-9; gotoxy(47,10); readln(rpm); gotoxy(47,11); readln(boostpressure); gotoxy(47,12); readln(deltapressure); gotoxy(47,13); readln(airtemp); gotoxy(47,14); readln(load);l triggercount:=O; rpm:=OoO; boostpressure:=OoO; deltapressure:=OoO; airtemp:=OoO; load:=OoO; gotoxy(47,9); readln(directionmeasure); pi:=4oO*atan(1o0); fringespacing:=488oOE-9/(2oO*sin(angle*pi/400o0)); lpossible velocitiesl preshiftvelocity:=2oO*numoflines*freqofgrating•fringespacing; const1:=numofchan-2; const2:=(fringespacing/tau)/(2o0*numofchan•(numofchan-1)); for i:=firstindex to velnumofchan do
possvelocity[i]:=(const1*(i+firstchan-1)+numofchan)*const2-preshiftvelocity;
numofexp:=O
1---l lthis procedure controls the correlatorl
procedure intakecorr(var autocorrelation :corvector; :integer; :integer; :integer; :integer; :integer}; var firstchan firstindex var intfail lengthcorr var orderotmatrix label 99; const var reset =1; start =2; seqmode =3; stop =4; command continue correlator defaultcammand i validcouand :char; :boolean; :rawdata; :char; :integer; :boolean;
lactual correlator output}
begin lof procedure intakecorrl continue:=true;
defaultcommand:='l'; while continue do begin
repeat (loop to get a valid commandl gotoxy(O,Ol;
write('lntake command: R(eset, S(tart, sT(op, I(ntake, ' 'A(uto, Q(uit; default=',defaultcommand,' '); gotoxy(70,0l;
readln(command);
if (ord(command)=32) or (command='R'l or (command='S'l or (command='T'l or (command='l'l or· (command='Q'l or (command='A'l then validcommand:=true el se validcommand:=false until validcommand; if ord(command)=32 then command:=defaultcommand; case command of 'R' :begin OPERATEK7023(reset); defaultcommand:='S' end; 'S' :begin OPERATEK7023(start); defaultcommand:='T' end; 'T' :begin OPERATER7023(stop); defaultcommand:='l' end; 'I' :begin end;
writeln('press any key to stop'); READK7023(correlator,lengthcorr,intfaill; for i:=firstindex to orderofmatrix-1 do
autocorrelation[i]:=correlator[i+3+firstchan]; autocorrelation[orderofmatrix]:=O.O; if intfail<>O then goto 99 el se defaultcommand:='Q' 'A' :begin end; OPERATER7023(reset); OPERATEK7023(start);
writeln('press any key to stop'); READK7023(correlator,lengthcorr,intfaill; for i:=firstindex to orderofmatrix-1 do
autocorrelation[i]:=correlator[i+3+firstchan]; autocorrelation[orderofmatrix]:=O.O; if intfail<>O then goto 99 el se defaultcommand:='Q' 'Q' :continue:=false
end lof case statement! end; lof while continue loop} 99: exit(intakecorrl
begin
writeln('enter meanvelocity'); readln(p1);
writeln('enter standarddev'); readln(p2);
writeln('v perpendicular to meanvel. '); readln(p3);
writeln('enter visibility {\)'); readln(p4);
writeln('enter noiselevel (\)'); readln(p5);
writeln('enter number of fringes 1 ·~p6:2:11 'in matrix');
readln(p6) end; {of procedure)
{---)
{this procedure generates a theoretica! autocorrelationl procedure generatecorr(var autocorrelation :corvector;
var arg1 arg2 exp1 exp2 j realrand
var firstchan :integer;
var fringespacing :real;
var numoffringes :real;
var orderotmatrix :integer;
var tau :real;
var meanvelocity :real;
var standarddev :real;
var perpendvel :real;
var visibility :integer;
var noiselevel :integer);
:real; :real; :real; :real; :integer; :real; begin for j:=1 to orderofmatrix-1 do begin arg1:=-(meanvelocity•meanvelocity+perpendvel•perpendvel)• (j+firstchan-1)*(j+firstchan-1)*tau•tau/ (fringespacing*fringespacing•numoffringes•numoffringes); i f arg1>-50 then exp1:=exp(arg1) else exp1 :=0.0; arg2:=-standarddev•standarddev•(j+firstchan-1)•(j+firstchan-1)• tau•tau*19.73921/(fringespacing•fringespacing); i f arg2>-50 then exp2:=exp(arg2) el se end fringespacing))+20.0)*1E6; realrand:=(random-16384.0)/32767.0; autocorr~lation[j):=autocorrelation[j]+ (1.0+visibility/100.0)•realrand•noiselevel•1E4; write( •. 1 )
var i ma x min xplot yplot :integer; :real; :real; :integer; :integer;
var firstchan :integer;
firstindex :integer;
var orderotmatrix :integer);
begin lof procedure sketchautol max:=autocorr[2];
min:=autocorr[2];
for i:=2 to orderofmatrix-1 do begin if autocorrelation[i]>max then max:=autocorrelation[i]; if autocorrelation[i]<min then min:=autocorrelation[i] end; page(output); for i:=1 to 72 do begin end; gotoxy(i+6,21); write( '-') for i:=O to 6 do begin end; gotoxy(12*i+6,22); write(12*i+firstchan:2) gotoxy(38,23l; write('chan. number'); for i:=O to 21 do begin gotoxy(6,i); write('! 'l end; for i:=firstindex to 72 do beqin xplot:=i+7; yplot:=trunc(21-20*(autocorrelation[i]-min)/(max-min)); gotoxy(xplot,yplot); write( '+' l end; i:=O; while min > 10 do beqin min:=min/10.0; i:=i+1 i:=O; while max > 10 do begin end; max:=max/10.0; i:=i+1 gotoxy(O, 1 l; write(max:1:1, 'E',i) end; lof procedure sketchautol
label 99; con st lowbound upbound var background beamradius const1 const2 const3 i j ma x mean numofpoints pi taufreq
firstindex :integer; (parameters)
var freqofgrating :real;
var fringespacing :real;
var invmatrix :matrix;
var jmeansq :real;
:integer; :integer; :real; :integer; :vel vector; :boolean; :real; :real; :integer; var moddepth var numofchan var numoffringes var orderafmatrix var passvelocity var scaled var tau var taumeansq velnumofchan
var intfail :integer; (output
var peakindex :integer; (parameters)
var velocprob :vel vector l ;
=5; =65; :real; :real; :real; :real; :real; :integer; :integer; :real; :real; :integer; :real; :vel vector;
lpeak search of velocity distributionl lin the interval lowbound to upboundl
(dummy constants to speed ( up the calculations (loop index)
(loop index!
lthe maximum value of the velocity probability vector) (mean value of autocorrelationl (number of points to stroke
the exponential fit! 13.1417 ... )
(vector of possible frequencies multiplied by tau)
begin lof procedure multiplycorrl
beamradius:=numoffringes*fringespacing; pi:=4.0*atan(1.0);
leerreetion if a scaled matrix is usedl
i f scaled then begin
const1:=taumeansq/(numoffringes•numoffringes); for j:=firstindex to orderofmatrix-1 do
{multiplication of autocorrelation with inverted matrix!
mean:=O.O; i:~O;
for j:=firstindex to orderofmatrix-1 do mean:=mean+autocorrelation[j]; mean:=mean/(orderofmatrix-1);
while mean > J.OE4 do
begin
mean:=mean/10.0; i:=i+1 end;
mean:=trunc(mean)*pwroften(i);
for j:=firstindex to orderofmatrix-1 do
autocorrelation[j]:=autocorrelation[j]-mean; writeln('start of multiplication');
(*for i:=firstindex to velnumofchan do begin
velocprob[i]:=O.O;
for j:=firstindex to orderafmatrix do velocprob[i]:=velocprob[i]+
invmatrix[i,j]•autocorrelation[j]; write('.' l
end; {of for i-statement)*)
matrixmultiply(invmatrix,autocorrelation,velprob, velnumofchan,orderofmatrix); for j:=firstindex to orderofmatrix-1 do
autocorrelation[j):=autocorrelation[j]+mean; i f scaled then begin const1:=numofchan-2; const2:=2.0*numofchan•(numofchan-1l; const3:=jmeansq/(numoffringes•numoffringes); for i:=firstindex to velnumofchan do
begin
taufreq[i]:=(const1*(i+1)+numofchan)/const2; velocprob[i]:=velocprob[i]•
exp(taufreq[i]•taufreq[i]•const3) end lof for statement!
end; lof if statement!
(calculation of background, only possible if the last matrixrow is readl lbackground:=O.O;
for j:=firstindex to orderafmatrix do
background: =background+ .
invmatrix[orderofmatrix,j]•autocorrelation[J];I lnormalise the peak of the velocity distributionl
writeln;writeln;
writeln('peak search and normalisation'); peakindex:=lowbound;
max:=0.01;
begin
if velocprob[i-1]>-0.0S*velocprob[i] then {this is supposed to rule out oscillation peaksl begin
if velocprob[i]>max then begin
max:=velocprob[i]; peakindex:=i end {of inner ifl end lof outer ifl end; lof for i statement) if max>0.01 then
el se
for i:=firstindex to velnumofchan do velocprob[i):=velocprob[i]/max begin
intfail:=300;
goto 99 {transfer to end of procedure) end; lof else clausel
99: exit(multiplycorr)
end; lof procedure multiplycorrl
1---l lthis procedure sketches the velocity probability distribution.l
procedure sketchvel{ firstindex :integer;
var
i
x plot yplot
var passvelocity :velvector;
var velocprob :velvector);
:integer; :integer; :integer;
{loop index)
begin lof procedure sketchveil page{output); gotoxy{5,21); write('---0 r '---'); gotoxy(2,22); write(possvelocity[1]:2:1); for i:=1 to 6 do begin gotoxy(12*i+2,22); write(possvelocity[12*i-1]:2:1) end; gotoxy(30,23); write('velocity (m/s)'); for i:=O to 21 do begin gotoxy( 4 r i); write('! ') end; for i:=O to 5 do begin gotoxy(O, 4*i+1); write((S-i)/5.0:1:1) end; gotoxy(0,7); write( 'P(v) '); for i:=firstindex to 72 do begin end
if (velocprob[i)>-1.0/20) and (velocprob[i]<22.0/20.0) then begin xplot:=i+S; yplot:=trunc(21-velocprob[i)*20); gotoxy(xplot,yplot); write( '+') end
1---l
lthis procedure calculates the meanvelocity, standarddev. and the squared differences between data and fit}
procedure expf i t ( f irs tindex :integer; var peakindex :integer;
var xdata :velvector;
var ydata :velvector;
var curtosis :real;
var meanvelocity :real; var numofpoints :integer;
var skewness :real;
var standarddev :real); const maxnum =8; dim =3; var i :integer; low :integer; moment1 :real; moment2 :real; moment3 :real; moment4 :real; pointsleft :integer; pointsright :integer; SUIIXY :real; sumx2y :real; sumx3y :real; sumx4y :real; sumy :real; x :real; y :real; begin lofmainl if peakindex-maxnum>firstindex+1 then low:=peakindex-maxnum el se low:=firstindex+1; i:=peakindex; repeat i:=i-1; pointsleft:=peakindex-i until (ydata[i)<0.1) or (i<=low); i:=peakindex;
repeat i:=i+1;
pointsright:=i-peakindex
until (ydata[i)<0.1) or (i>=peakindex+maxnum); numofpoints:=pointsleft+pointsright+1; sumy:=O.O; sumxy:=O.O; sumx2y:=O.O; tinput I (output) sumx3y:=O.O sumx4y:=O.O
for i:=peak ndex-pointsleft to peakindex+pointsright do begin end; x: =xdata[i); y:=ydata[i]; sumy:=sumy+y; sumxy:=sumxy+x*y; sumx2y:=sumx2y+x*x*y; sumx3y:=sumx3y+x*x*x*y; sumx4y:=sumx4y+x*x*x*x*y moment1:=sumxy/sumy; moment2:=sumx2y/sumy; moment3:=sumx3y/sumy; moment4:=sumx4y/sumy; meanvelocity:=moment1; standarddev:=sqrt(abs(moment2-moment1*moment1}); skewness:=(moment3-3.0*moment1*moment2+2.0*moment1*moment1*moment1)/ standarddev•standarddev*standarddev; curtosis:=(moment4-4.0*moment1*moment3+6.0*moment1*moment1*moment2-3.0*moment1*moment1*moment1*moment1)/ standarddev*standarddev*standarddev*standarddev end; lof procedure expfitl
end; end el se defaultnext:='I' 'F' :begin expfit(frstindex,ipeak,possvel,velprob,curtos, meanvel,numpoints,skewnes,standev); writetoscreen(curtos,numpoints,meanvel, skewnes,standev); defaultnext:='S' end; 'S' :begin writetofloppy(airtem,autocorr,boostpres,curtos,day, deltapres,dimeasure,frstchan,frstindex, freqgrating,frinspacing,ifail, lstindex,lload,meanvel,month,numchan, numexp,numpoints,rradius,rev, skewnes,standev,ttau,trigcount); end; writetoprinter(autocorr,boostpres,curtos,day, deltapres,dimeasure,lload, meanvel,month,numexp,rradius,rev, skewnes,standev,trigcount,velprob); lclose(printerfile); close(diskfile,lockl;l defaultnext:='I' 'C' :begin end;
page(output); lclear screen} gotoll:y(0,2l; contrparameter(frstchan,freqgrating,frinspacing,mdepth, numchan,numfringes,scale, ttau) 'G' :begin end; page(output); param(meanvel,standev,perpvel,visib,noise,numfringes); generatecorr(autocorr,frstchan,frinspacing,numfringes, ordermatrix,ttau,meanvel,standev, perpvel,visib,noise); defaultnext:='M' 'E' :begin end; page(output); readcorr(autocorr,fileid,frstindex,ifail,lstindex, numchan,ordermatrix,ttau); defaultnext:='M' 'T' :begin continue:=false; close(printerfile); close(diskfile,lock); 99: if ifail<>O then errorhandler(blcknumber,fileid,ifail)
en~; lof whi~e continue loop) ent(program)
