##################################################################
# SIMUL.WS Workspace zur Simulation der 3D Laser-Abklingkurve #
##################################################################
# Groesse des Simulationsvolumens (In Bilddefinitionen zu aendern)
# dx := 101; dy := 61; dz := 200;
# Simulationsvolumen
float simvol[200][61][101];
# Kamerabild
float image[61][101];
# ZT(X)-Bild (Zum Abspeichern der zeitlichen Entwicklung des Tiefenprofiles)
float zt[200][1][123];
# Temperaturerhoehung im Volumen durch Laser (3-D LUT)
float laser[200][61][101];
# Hilfsbilder zum Erzeugen der Laserintensitaet
float row[101]; # x-Ausdehnung des Simulationsvolumens 1 pixel = 1 mm
float col[61][1]; # y-Ausdehnung des Simulationsvolumens 1 pixel = 1 mm
float depth[200][1][1]; # z-Ausdehnung des Simulationsvolumens 1 pixel = 10 micron
# Bilder zur Berechnung der Volumen-Abstrahlung
# data: Daten als ascii File abgespeichert, muss eingelesen werden
# y:i, x: [0]lambda_i [1]B_i(288) [2]B_i(228) [3]deltalambda_i [4]zeta_i
# volflux: additive Temperaturaenderung in z-Richung (berechnet sich aus data)
float data[74][5];
float volflux[200][1][1];
# image-LUT zum Abspeichern der tiefenabhaengigen und Temperaturabhaengigen
# Werte, die zum Entstehen des Kamerabildes beitragen.
# z: Tiefe x:Temperatur 12,...,15,...31 Grad Celsius
float B_jk[200][1][20];
# Faltungsmasken (allg. [k,1-2k,k], in horizontaler Richtung ist k um den Faktor 10^4
# kleiner, da die Ortsaufloesung um den Faktor 100 groesser ist)
# SURFACE RENEWAL MODELL -> Die Summe der Koeffizienten ist 1 - dt/tau, da Waerme verschwindet,
# mit tau = 2 sec und dt = 1.78E-4 s -> dt/tau = 8.9E-5.
# NUR FUER DIE VERTIKALE FALTUNG!
binh := {0.000025,0.99995}; # horizontal
#binv := {0.25,0.5}; # vertikal (reine Diffusion)
binv := {0.25,0.499911}; # vertikal (surface renewal)
# LASER: Erzeugen der additiven Temperaturzunahme durch den Laser
# leistung = Gesamtleistung [Watt]
# sigmasquare = Varianz der horiz. Ausdehnung des Lasers (sigma^2) [mm^2]
# zeta = Eindringtiefe der Laserstrahlung [microns]
operator laser(leistung,sigmasquare,zeta);
# Berechnung der Normierungskonstante (2 Pi Sigma^2 zeta)^-1
float normierung;
normierung = leistung;
normierung = Div(normierung,6.3); # 2 * Pi
normierung = Div(normierung,sigmasquare);
normierung = Div(normierung,zeta);
# Erzeugen eines Zeilen und Spalten-Gauss
row = Wedge({-50,1.0});
col = Wedge({-30,1.0,0.0});
row = Sqr(row); row = Neg(row);
col = Sqr(col); col = Neg(col);
# Standardabweichung 2 sigma^2
row = Div(row,sigmasquare); col = Div(col,sigmasquare);
row = Div(row,2.0); col = Div(col,2.0);
row = Exp(row); col = Exp(col);
# 2-D Gauss (konstant in z) durch Multiplikation der Zeilen und Spalten
laser = row;
laser = Mul(laser,col);
# Erzeugen der exponentiellen Abklingkurve in z-Richtung
depth = Wedge({0.0,10.0,0.0,0.0}); # Stufe 10.0, da dz = 10 microns
depth = Div(depth,zeta); depth = Neg(depth);
depth = Exp(depth);
depth = Mul(depth,normierung);
# Multiplikation des 3-D Gauss mit der Tiefenfunktion
laser = Mul(laser,depth);
# Umrechnen von Leistung in Delta-Temperatur -> [laser] = 1 K
laser = Mul(laser,42.7);
endoperator;
# VOLFLUX: Volumenemission langwelliger Strahlung
# erzeugt einen Vektor in z-Richtung, der die Temperaturabnahme
# pro Volumenelement der Simulation pro Zeitschritt enthaelt [K].
# Er muss bei jedem Zeitschritt von der Temperatur SUBTRAHIERT werden
# ACHTUNG: Dieser Wert gilt nur fuer die Temperaturen Twasser=15°C und Tsky= -45°C
# [Saunders,70]. Fuer andere Temperaturen muss der Wert fuer B_Twasser-B_Tsky
# neu berechnet werden.
# depthlamb: Hilfsbild mit Termen der Summe, die noch lambda-abhaengig ist (i)
# templamb: Hilfsvektor
operator volflux();
float depthlamb[200][74][1];
float templamb[74][1];
# Lesen der vorab berechneten Werte (siehe Definition v. data)
data = Read("c:\user\hhaus\longrad.asc");
depthlamb = Wedge({0.0,10.0,0.0,0.0}); # Stufe 10.0, da dz = 10 microns
depthlamb = Mul(depthlamb,1.5); # Korrektur fuer Tiefenabsorption
depthlamb = Div(depthlamb,data[][4]);
depthlamb = Neg(depthlamb);
depthlamb = Exp(depthlamb);
depthlamb = Div(depthlamb,data[][4]); # / zeta
depthlamb = Mul(depthlamb,10.0); # * dz
depthlamb = Mul(depthlamb,data[][3]);
templamb = Sub(data[][1],data[][2]);
depthlamb = Mul(depthlamb,templamb);
volflux = Sum(depthlamb);
# Umrechnen von Fluss j_k auf deltaT_k [K]
volflux = Div(volflux,234192.0);
endoperator;
# RANDOBEN: Randbehandlung an der Wasseroberflaeche. Die Faltungsformel
# fuer den oberen Rand hat die Form C0 = 2k*C1 + (1-2k)*C0 + 2(dt/dz)(j/ro/c),
# die Faltung fuer den Bulk die Form C0 = k*C1 + (1-2k)*C0 + k*C-1. Da aber
# bei der Faltung c-1 = 0 angenommen wird, ist der Randterm nicht falsch, sondern
# nur zu klein um den Betrag k*C1 + 2(dt/dz)(j/ro/c). Dieser muss hinzuaddiert werden.
# fluss = Flussdichte an der Wasseroberflaeche [W/m^2] = j_evaporation + j_sensible.
operator randoben(fluss);
float fbild[61][101];
float addfluss;
addfluss = Mul(fluss,0.00854);
addfluss = Mul(addfluss,0.001);
fbild = Mul(simvol[1],0.25);
fbild = Sub(fbild,addfluss);
simvol[0] = Add(simvol[0],fbild);
endoperator;
# RANDUNTEN: Randbehandlung an der Grenze zum bulk.
# Die Bulktemperatur wird als konstant angenommen. Daher wird der unterste
# Kasten (Bild) nach jedem Simulationsschritt, d.h. jeder Faltung in z-Richtung
# wieder auf die Temperatur T0 gesetzt. Dadurch ragt die Faltungsmaske fuer
# das vorletze Bild immer in die Temperatur T0.
operator randunten(T0);
simvol[199] = T0;
endoperator;
# RANDAUSSEN: Randbehandlung an der vertikalen Berandung des Simulationsvolumens.
# Es werden zyklische Randbedingungen angenommen. Die Temperaturen an den Aussen-
# grenzen des Simulationsvolumens werden daher denen an den gegenueberliegenden
# Waenden gleichgesetzt. Und zwar immer dem vorletzten vor der Wand, da die aeusseren
# durch die Faltung mit 0 verfaelscht werden.
operator randaussen();
simvol[][][0] = simvol[][][99];
simvol[][][100] = simvol[][][1];
simvol[][0][] = simvol[][59][];
simvol[][60][] = simvol[][1][];
endoperator;
# IMAGELUT: Berechnen der Verschiebe-LUT zur Berechnung der Bildintensitaet
# Das Bild B_jk enthaelt in z-Richtung den Beitrag dieser Tiefenschicht zum
# Kamerabild. Da nicht a priori bekannt ist, welche Temperatur in dieser Tiefe
# vorliegt, ist die LUT temperaturabhaengig (in x-Richtung). Um das Kamerabild
# an einem Punkt zu erhalten, muss die LUT:
# 1. verschoben werden, so dass die richtigen Temperaturen uebereinander liegen und
# 2. ueber die Tiefe (in z-Richtung) aufaddiert werden
# Die dazu noetige Verschiebung betraegt T(z)-T0. Danach stehen die richtigen
# Temperaturen in B_jk[][3]
operator ImageLUT();
float camera[10][22];
float B_ij[10][20];
float B_ijk[200][10][20];
float e_jk[200][10][1];
# Lesen der vorab berechneten Werte fuer zeta_j, s_j, und B(j,T_i)
camera = Read("c:\user\hhaus\camera.asc");
# Umkopieren der Werte B(j,T_i) in ein 3D-Bild
B_ij = camera[][2:]; B_ijk = B_ij;
# Erzeugen des tiefenabhaengigen Absorptionskoeffizienten, wellenlaengenabhaengig
e_jk = Wedge({0.0,10.0,0.0,0.0}); # Schrittweite 10, da dz = 10 microns
e_jk = Div(e_jk,camera[][0]); # camera[][0] enthaelt die Werte zeta_j
e_jk = Neg(e_jk); e_jk = Exp(e_jk);
e_jk = Div(e_jk,camera[][0]); # Faktor 1/zeta
# Multiplikation der Extinktionswerte mit den spektralen Empfindlichkeiten
e_jk = Mul(e_jk,camera[][1]);
# Multiplikation der Planckkoeffizienten mit der Tiefenabhaengkeit und Empfindlichkeit
B_ijk = Mul(B_ijk,e_jk);
# Integration (Summe) ueber die Wellenlaengen
B_jk = Sum(B_ijk);
endoperator;
# IMAGE: Erzeugt aus der gegebenen Temperaturverteilung im 3D Simulationsvolumen
# und der tiefen- und temperaturabhaengigen Planckterme die Intensitaets-
# verteilung auf dem Kamerabild (image).
# B_jk muss vorher mit ImageLUT() erzeugt werden.
operator Image();
float shift[200][1][1];
float temp[200][1][101];
float B_jk_temp[200][1][20];
image = 0.0;
scan(image,y,simvol,y);
temp = simvol;
scan(image,x,temp,x);
shift = temp;
shift = Sub(temp,15.0); # Temperaturdifferenz um dT=T(z)-T_0
B_jk_temp = Shift(B_jk,shift); # Verschieben der Temperaturen um dT
image = Sum(B_jk_temp[][0][3]); # Aufaddieren der Tiefenbeitraege an
endscan; # der Stelle T_0+dT
endscan;
endoperator;
# INIT_SIM: Hier werden die LUTs berechnet und die Startwerte initialisiert
operator init_sim();
simvol = 15.0;
laser(17.3,25.0,11.5);
showlaserxz();
volflux();
ImageLUT();
endoperator;
# SHIFT: Schiebt den Inhalt des Simulationsvolumens um V_SURF [Pixel] in positive
# x-Richtung. Nur positive shifts bis zu 101, zyklisch, d.h. ein shift
# um 101 ergibt das Originalbild, wenn man einen negativen shift der
# Groesse dx will, dann muss man um 101-dx nach rechts schieben.
# Die Groesse DV_DZ gibt den Gradienten des Geschwindigkeitsfeldes an.
# Pro Tiefenschritt (10 micron) wird die Verschiebung um DV_DZ [pixel/Schritt]
# kleiner (VORSICHT: Bei zu grossem Gradienten kann die Verschiebung ab einer
# bestimmten Tiefe negativ werden)
operator shift(v_surf,dv_dz);
float tempimg[61][202];
float shift[200][1][1];
float gradient[200][1][1];
float wz[4];
wz = 0.0; wz[1] = dv_dz;
shift = 101.0;
gradient = v_surf;
shift = Sub(shift,gradient);
gradient = Wedge(wz);
shift = Add(shift,gradient);
scan(simvol,z,shift,z);
tempimg = simvol; tempimg[][101:] = simvol;
simvol = Shift(tempimg,shift);
endscan;
endoperator;
# SIMULATION: Dies ist der Operator, der die Simulation durchfuehrt
operator simulation();
short bildnummer;
string snummer;
string filename;
bildnummer = 0;
SetFormat(bildnummer,"%03d");
Image();
Write("c:\user\hhaus\nullbild.raw",image);
repeat(3); # 3 * 16.66 ms
repeat(93); # LASER AN, 93 Simulationsschritte = 16.66 ms = 1 Halbbild
simvol = Add(simvol,laser); # Temperaturerhoehung durch Laser
simvol = Sub(simvol,volflux); # Temperaturerniedrigung durch longwave emission
shift(0.017857,0.000105); # Konvektion mit Scherstroemung
simvol = ConvSep(simvol,binv,"zs"); # Diffusion vertikal
simvol = ConvSep(simvol,binh,"xsys"); # Diffusion horizontal
simvol = Add(simvol,0.001335); # Korrektur ... + dt/tau*C_0 fuer Surface Renewal Modell
randoben(150.0); # Randbedingung oben, Oberflaechenfluss 150 W/m^2
randunten(15.0); # Randbehandlung unten, bulk Temperatur
randaussen(); # periodische Randbedingungen
endrepeat;
zt[][0][bildnummer] = simvol[][30][50]; # Tiefenschnitt in zt-bild speichern
Image();
bildnummer = Add(bildnummer,1);
endrepeat;
repeat(10); # 12 * 10 * 16.66 ms = 2 s
repeat(12); # Alle 12 Halbbilder = 0.2 s wird das Simulationsvolumen abgespeichert
repeat(93); # LASER AUS, 93 Simulationsschritte = 16.66 ms = 1 Halbbild
simvol = Sub(simvol,volflux); # Temperaturerniedrigung durch longwave emission
shift(0.0,0.000105); # Konvektion mit Scherstroemung
simvol = ConvSep(simvol,binv,"zs"); # Diffusion vertikal
simvol = ConvSep(simvol,binh,"xsys"); # Diffusion horizontal
simvol = Add(simvol,0.001335); # Korrektur ... + dt/tau*C_0 fuer Surface Renewal Modell
randoben(150.0); # Randbedingung oben, Oberflaechenfluss 150 W/m^2
randunten(15.0); # Randbehandlung unten, bulk Temperatur
randaussen(); # periodische Randbedingungen
endrepeat;
zt[][0][bildnummer] = simvol[][30][50]; # Tiefenschnitt in zt-bild speichern
Image();
bildnummer = Add(bildnummer,1);
endrepeat;
endrepeat;
endoperator;