Numerik: Wegintegral lösen

Jester

agga ugga schrieb:

naja, das e^(-T) scheint ja eine konstante zu sein. wenn nicht, dann ist die formulierung des integranden ein wenig ungeschickt.
wenn e^(-T) konstant ist, kann man es vor das integral ziehen und es gibt eine einfache stammfunktion -e^(-h/8000)

Lass sehen.
Bitte denk dran, dass da ein h steht, Du aber entlang einem Weg r läufst, der nicht senkrecht nach oben geht.

this-*gt*that

Hallo Leute

da mein Problem doch eher zu verwirren scheint, werd ich es nochmal verkürzt schildern.

Ich will die sogenannte Optical Depth (T) berechnen (ein Maß, dass die "Lichtundurchdringlichkeit eines Mediums wie Luft" beschreibt. T << 1 heißt optisch dünn und T >> 1 optisch dick (z.B. eine Wolke oder verschmutzter Himmel )

Diese Optical Depth T lässt sich berechnen per:
T = \frac{4\*Pi\*Kr}{L^4}\int_{P\_a}^{P\_b}e^\frac{-h}{H_0} ds
Der Bruch vorm Integral ist im Grunde nur eine Konstante, was ich berechnen will ist nur der Integralteil. Dieses $e^\frac{-h}{H_0}$ beschreibt die Dichte einer Höhe h über Erdboden. H_0 ist eine Konstante (z.B. 8000m).
Als Koordinatensystem hab ich nach wie vor das: http://www.infoboard.org/screenshots/system.png
Zusammengefasst: Ich muss also die Dichte entlang der Strecke zwischen den Punkten P_a und P_b (der blaue Strich) integrieren.

Dazu habe ich diesen Sourcecode gefunden:

double height = EPSILON + planetaryRadius + x; // x ist zwischen 0 und 80 000
D3DXVECTOR3 pos(0, height, 0);  // pos ist am Anfang P_a

double distance = ... // distance = Länge der gesamten Strecke (blau im Bild)

double sampleLength = distance / samples; // Länge eines Sample Intervalls (in realen Metern)

double scaledLength = sampleLength * scale;  // Normalisierte Länge eines Sample Intervalls. scale = 1 / 80000. 

D3DXVECTOR3 sampleRay = ray * sampleLength;
pos += sampleRay * 0.5f;  // pos ist der jeweilige Sample Point entlang der Strecke (blau)		       
double opticalDepth_R = 0;
// Über alle Sample Points iterieren um die Optical Depth aufzusummieren 
for(int i=0; i< samples; i++) {
   double height = D3DXVec3Length(&pos); 
   double samplePointAltitude = (height - planetaryRadius) * scale; 

   opticalDepth_R += exp(-samplePointAltitude / 0.1);
   pos += sampleRay;
}
opticalDepth_R *= scaledLength;

Alles was der Code macht ist also über ein paar Sample Points entlang der Strecke gehen, e^(-h/H0) auswerten und aufsummieren und dann am Ende mit scaledLength multiplizieren.

Aber ich verstehe nicht, wieso er mit scaledLength multipliziert.
Erst dachte ich, er stellt einfach Rechtecke unters Integral, aber die Breite eines Recktecks wäre dann ja sampleLength. Und eigentlich kann er ja nicht einfach so Rechtecke drunterstellen, weil es ja kein normales Integral sondern ein Wegintegral ist. Also müsste dieses scaledlength eher die 1. Ableitung der Parametrisierung sein, aber das ergibt auch keinen Sinn...
Wieso wird also mit scaledLength multipliziert?

Daniel E.

this->that schrieb:

Erst dachte ich, er stellt einfach Rechtecke unters Integral, aber die Breite eines Recktecks wäre dann ja sampleLength. Und eigentlich kann er ja nicht einfach so Rechtecke drunterstellen, weil es ja kein normales Integral sondern ein Wegintegral ist. Also müsste dieses scaledlength eher die 1. Ableitung der Parametrisierung sein, aber das ergibt auch keinen Sinn...
Wieso wird also mit scaledLength multipliziert?

Also flapsig gesprochen passiert doch folgendes:
Für ein Wegintegral nimmst Du den Funktionswert an jeder Stelle und summierst die zusammen. Wenn man jetzt so tut, als wäre die Funktion über eine kurze Strecke ds konstant, kann man den Funktionswert an einer Stelle mit der Länge der Strecke ds gewichten und ds ist nunmal |r'| dt, mit r als deiner Parametrisierung.

Die Unterscheidung "normales Integral" und Kurvenintegral ist keine so schrecklich fundamentale, sobald man eine Parametrisierung festgelegt hat, kannst Du natürlich auch "Rechtecke drunterlegen".

Welche Parametrisierung hat der Autor also festgelegt? Soweit ich das sehe, läuft er den Strahl einfach von A nach B, also ist die natürliche Parametrisierung die Bogenlänge, also ist |r'|=1, bzw. ds=dt, und jetzt approximiert man normalerweise die Sache durch eine Summe, so daß ds~=Abstand/N.

Woher der Skalierungsfaktor kommt, ist also eine berechtigte Frage und gar nicht klar. Ich vermute, daß es nichts mit der Integration an sich zu tun hat, sondern irgendwelche Einheiten oder Faktoren umnormiert. Könnte das evtl sein? Ist "deine" Integration genau um den Faktor daneben?

Jester schrieb:

Lass sehen.
Bitte denk dran, dass da ein h steht, Du aber entlang einem Weg r läufst, der nicht senkrecht nach oben geht.

naja, wenn ich mich nicht irre, dann braucht man bloss den vektor pa auf pb zu projizieren und dann einfach vom projizierten pa bis zum betrag des vektors pb zu integrieren. damit läuft man quasi entlang einer senkrechten nach oben.

rüdiger

this->that schrieb:

Der Sourcecode, den ich gefunden habe, sieht so aus:
float scale = 1 / 80000; // 1 / Atmosphärenhöhe. Normalisiert Höhen in [0,1] 
float sampleLength = distance / sampleCount; // Länge eines SampleIntervalls (distance = StreckenLaenge, also Länge der blauen Strecke)
float scaledLength = sampleLength * scale;   // Skalierte SampleLänge, eine Länge von 80000 wäre also 1
float3 sampleRay = ray * sampleLength;
float3 samplePos = camPos + sampleRay * 0.5f;  // jeweilige Sample Position entlang des blauen Strahls
    
float sum = 0;

for(int i=0; i < sampleCount; i++) { 
   float h = length(samplePos) - innerRadius;	
   sum += exp(-h/8000)*exp(-T);
		
   samplePos += sampleRay;	
}
float endResult = scaledLength * sum;
Im Grunde macht der Code also nichts anderes als: Iteration über ein paar SamplePoints entlang der Strecke, bei jedem SamplePoint die Funktion auswerten und zum nächsten SamplePoint gehen.
Ganz am Ende wird dann mit der skalierten Länge eines SampleIntervalls multipliziert.

Diesen Code verstehe ich nicht, aber er funktioniert.
Er geht ja ganz normal die echte Strecke entlang, wieso multipliziert er dann am Ende mit scaledLength? Benutzt der überhaupt ein Wegintegral?

Nicht das ich mich gut in dem Bereich auskenne. Aber das schaut mir einfach nach einem Riemannintegral aus.

Man wertet die Funktion an den Samplepositionen aus, die man entlang des Strahls verschiebt und multipliziert sie mit der Breite. Da man die Verschiebung entlang des Strahls macht, kann man ja durch Rechtecke nähren.

(Disclaimer: Ich lerne das Thema gerade erst selbst und das Integral-Kapitel kommt erst die nächsten Tage dran.)

Jester

agga ugga schrieb:

Jester schrieb:

Lass sehen.
Bitte denk dran, dass da ein h steht, Du aber entlang einem Weg r läufst, der nicht senkrecht nach oben geht.

naja, wenn ich mich nicht irre, dann braucht man bloss den vektor pa auf pb zu projizieren und dann einfach vom projizierten pa bis zum betrag des vektors pb zu integrieren. damit läuft man quasi entlang einer senkrechten nach oben.

ne, damit veränderst du den weg den du läufst... und das darfst du nur wenn es eine stammfunktion gibt.

Jester schrieb:

ne, damit veränderst du den weg den du läufst... und das darfst du nur wenn es eine stammfunktion gibt.

wenn das integral wegunabhängig ist, dann darf man das ja.
dann kommt es in diesem fall nur auf das Δh an.
dieses e^(-T) scheint ja ne konstante zu sein, dann ist ja alles in buddäär

Jester

agga ugga schrieb:

Jester schrieb:

ne, damit veränderst du den weg den du läufst... und das darfst du nur wenn es eine stammfunktion gibt.

wenn das integral wegunabhängig ist, dann darf man das ja.
dann kommt es in diesem fall nur auf das Δh an.
dieses e^(-T) scheint ja ne konstante zu sein, dann ist ja alles in buddäär

Deine Argumentation geht grad ein bißchen im Kreis: "Wenn ich die Projektion machen darf, dann krieg ich ne einfache Stammfunktion und ich darf die Projektion wirklich machen, weil's dann wegunabhäng ist, also hat's ne stammfunktion und ist wegunabhängig"

Ist es aber nicht wegunabhängig, dann darfst Du nicht einfach projizieren und bekommst keine einefache Stammfunktion... usw.

Ich denke nicht, dass das Integral wegunabhängig ist. Schreib Dir die Formel für die Höhe doch mal explizit auf, das ist irgendein häßlicher Term mit ner Wurzel drin.

Jester

rüdiger schrieb:

Nicht das ich mich gut in dem Bereich auskenne. Aber das schaut mir einfach nach einem Riemannintegral aus.

Natürlich. Wegintegral werden doch darüber definiert.

rüdiger

Jester schrieb:

rüdiger schrieb:

Nicht das ich mich gut in dem Bereich auskenne. Aber das schaut mir einfach nach einem Riemannintegral aus.

Natürlich. Wegintegral werden doch darüber definiert.

Riemannintegral stimmt ja auch nicht ganz, da man ja die Unterteilung nicht unendlich fein macht und auch nicht Ober-/Untersumme.

Aber Daniel E. hatte es ja eh schon erklärt. Die Skalierung dürfte doch dadurch zustande kommen, dass einfach skaliert gerechnet wird: double samplePointAltitude = (height - planetaryRadius) * scale;

Jester

rüdiger schrieb:

Jester schrieb:

rüdiger schrieb:

Nicht das ich mich gut in dem Bereich auskenne. Aber das schaut mir einfach nach einem Riemannintegral aus.

Natürlich. Wegintegral werden doch darüber definiert.

Riemannintegral stimmt ja auch nicht ganz, da man ja die Unterteilung nicht unendlich fein macht und auch nicht Ober-/Untersumme.

Das tut der Richtigkeit der Aussage, dass Weg-Integrale über Riemann-Integrale definiert werden keinen Abbruch.

Aber Daniel E. hatte es ja eh schon erklärt. Die Skalierung dürfte doch dadurch zustande kommen, dass einfach skaliert gerechnet wird: double samplePointAltitude = (height - planetaryRadius) * scale;

Das sehe ich nicht. Das passiert zwar im Exponenten, warum es aber auch als Skalierfungsfaktor nochmal auftritt ist mir nicht klar.

ich habe nachgesehen, wie man mathematisch nachweist, das ein wegintegral vom weg unabhängig ist.
dafür muss man partielle ableitungen bilden, wenn die gleich sind, ist das wegintegral unabhängig.
kurzum:
nach wenigen zeilen wird klar, das dieses integral nicht vom weg unabhängig ist. :p