SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
     1 JLTYP,TEMP,PRESS,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION COORDS(3),FLUX(2),TIME(2)
      CHARACTER*80 SNAME
C ###############################################################
C Definition allgemeine Parameter
C ###############################################################

	PI=3.1415

C Oberflächenwärmequelle
	JLTYP=1

C Einlesen der aktuellen Koordinaten aus der Simulation
	x=COORDS(1)
	y=COORDS(2)
	z=COORDS(3)

C Einlesen der seit Simulationsstart verstrichenen Zeit
	dt=TIME(1)

C ###############################################################
C Definition Bauteileigenschaften
C ###############################################################

C Lage der Trennebene => Dicke des durchstrahlten Bauteils
C	z_Trenn=0.0005	!	z_Trenn=0.0005

C Absorptionskoeffizienten für die Bauteile und Inhomogenitäten
C Angabe in "1/m" 
C alpha = 8 entspricht alpha von PC Natur für Nd:YAG
C alpha = 10000 entspricht alpha von 0,1 Gew-% Ruß
	alpha_oben=20
	alpha_unten=20
	alpha_inhom=20
	alpha_lokal=20


C ###############################################################
C Definition Strahlparameter und Strahlkaustik
C ###############################################################

C Leistung des Strahls (in W)
	P=128

C Leistungsfaktor
	PGauss=2.39

C Formfaktor Gauss
	FG=2

C Strahlradius im Fokus senkrecht zur Strahlachse (in m)
	w_0=0.25E-3 !0.03035E-3

C Aperturwinkel Optik
	!zR=0.326E-3
	alphaAP=10.0!atan(sqrt(2.0)*w_0/zR)*180/PI
	BM_alphaAP=alphaAP*PI/180.0
C	write(6,*) alphaAP,BM_alphaAP
	
C Anfangskoordinaten der Fokusposition (in m)
	FP_x=0.000001
	FP_y=1.000E-3
	FP_z=4.000E-3


C Berechnung von Hilfsgrößen der Strahlgeometrie:
C Absolutwert des Abstands der aktuell berechneten Position vom Fokus entlang
C Strahlachse (in m)
	DIST_FP=ABS((FP_z-z))

C Strahlradius an Oberfläche
	w_0b=w_0+FP_z*TAN(BM_alphaAP)

C Strahlradius an aktueller z-Position senkrecht zur Strahlachse (in m)
	w_z=w_0+DIST_FP*TAN(BM_alphaAP)

C Flächenleistungsdichtemaximum an der Oberfläche (Faktor aus Gauß: 2,39 =
C (maximale Leistung) / (mittlerer Leistung» (in W/m²)
	A_0b=PI*w_z**2
	DP_0b=PGauss*(P/A_0b)

C ###############################################################
C Bewegung des Strahls mit Aktualisierung der Position
C ###############################################################

C Vorschubgeschwindigkeit (in m/s): ausschließlich
C x-Richtung
	v=0.5E-3
C y-Positionen konstant, da Bewegung ausschließlich in x-Richtung
C Berechnen der x-Positionen aus Anfangswerten, Geschwindigkeit und
C verstrichener Zeit
C	AFP_x=FP_x+v*dt
C	AFP_y=FP_y+v*dt

	FP_x=FP_x+v*dt
C	FP_y=FP_y+v*dt

C Berechnen des Abstandes der aktuellen Koordinaten vom Zentral strahl (in der
C x-y-Ebene)
C	dx=ABS(AFP_x-x)
C	dy=ABS(AFP_y-y)

	dx=ABS(FP_x-x)
	dy=ABS(FP_y-y)

C Berechnen des radialen Abstandes zum Zentral strahl (in der x-y-Ebene)
	r=SQRT(dx**2+dy**2)

C ###############################################################
C FLBB - FAKTOR LAMBERT-BEER-Gesetz
C ###############################################################

C Berechnen von Strahl
C an Einzelelement definierte Absorptionsänderung (Abbruchkriterium Exponent
C > 6)
C => Störstelle (exemplarisch Element Nr. 5081 im oberen Bauteil)

C ***********************************************
C AB HIER AUSKOMMENTIERT 
C ***********************************************
C	if (NOEL.eq.5081) THEN
C		ELBB=alpha_oben*z 
C		if (ELBB.le.6) THEN
C			FLBB=1/EXP(ELBB)
C		else
C			FLBB=0
C		endif
C			alpha_lokal=alpha_oben
C			alpha_lokal=alpha_inhom


C Absorption in Bauteil 1 (Abbruchkriterium Exponent > 6)
C	else if (NOEL.le.6000) THEN
	ELBB=alpha_oben*z
C	if (ELBB.le.6) THEN
C		FLBB=1.0/EXP(ELBB)
C	else
C		FLBB=0
C	endif
C	alpha_lokal=alpha_oben

C Absorption in Bauteil 2 (Abbruchkriterium Exponent > 6)
C	else
C		z_unten=ABS(z-z_Trenn)
C		ELBB=alpha_unten*z_unten+alpha_oben*z_Trenn
		if (ELBB.le.6) THEN
			FLBB=1/EXP(ELBB)
		else
			FLBB=0
		endif
C		alpha_lokal_=alpha_unten
C	endif

C ###############################################################
C FGAU - FAKTOR GAUSS-Verteilung
C ###############################################################

C Abbruch bei 3sigma; => 99,7% Abdeckung
	EGAU=FG*(r/w_z)**2
	if (EGAU.le.6) THEN
		FGAU=1/EXP(EGAU)
	else
		FGAU=0
	endif


C ###############################################################
C Berechnen der absorbierten Leistung an den aktuellen Koordinaten
C ###############################################################
C lokal auftreffende Flächenenergie [W/m']
	DPlokal=DP_0b*FLBB*FGAU

C bei Koordinaten absorbierte Leistung aus beiden Strahlen [W/m']
C => Übertrag in das Hauptprogramm
	FLUX(1)=DPlokal*alpha_lokal

C ###############################################################
C Skript-Ende
C ###############################################################
	
	
	RETURN
 	END
C ###############################################################
	
	 	SUBROUTINE FILM(H,SINK,TEMP,KSTEP,KINC,TIME,NOEL,NPT,
     1 COORDS,JLTYP,FIELD,NFIELD,SNAME,NODE,AREA)
C
	INCLUDE 'ABA_PARAM.INC'
C
	DIMENSION H(2),TIME(2),COORDS(3),FIELD(NFIELD)
	CHARACTER*80 SNAME

C—--Durchmesser und Wärmeübergangskoeffizient des gaußförmigen Luftstroms---
	Wo=3*4.0E-1 !war 0.0040E+3 ! "Taillenradius" gaußförm. Wärmeübergangsk.[mm]
	Wo2=3*1.6E-1 !war 0.0016E+3 ! "Taillenradius" vom abgezogenen
	! gaußförmigen Wärmeübergangskoeffizienten [mm]
	z_Mittelpunktwaermekoeffizient = 1800E-3 ! Maximum vom gaußförmigen
	! Wärmeübergangskoeff. [mW mm-2K-1]
	z_Mittelpunktwaermekoeffizient2 = 1300E-3 ! Maximum vom abgezogenen
	! gaußf. Wärmeübergangsk. [mW mm-2K-1] war 1300
	Streckungsfaktor = 1.8 ! Streckt das Gaußprofil ab dem
	! Mittelpunkt in x-Richtung

C----Definition vom Startpunkt und Endpunkt des Lufstrommittelpunkts und----
C----der 1. ,2. und 3. Streckenlänge bzw. der jeweiligen Geschwindigkeit----
C---(v2, v3, x2 und x3 wird nur bei variierender Geschwindigkeit benötigt)--
	xS = 0.0E+3 ! x-Startmittelpunkt vom Mittelpunkt des
	! gaußförmigen Wärmeübergangskoeffizient in [mm]
	
	xEnde = 0.2 ! Endpunkt vom Mittelpunkt des gaußförmigen
	! Wärmeübergangskoeffizient (bei Bewegung) in [mm]
	
	Abstand = 3*1.6E-3 ! x-Abstand der beiden Gaußmittelpunkte in [mm]

	xS2 = xS + Abstand ! x-Startpunkt vom Mittelpunkt des abgezogenen
	! gaußförmigen Wärmeübergangskoeffizienten in [mm]

	yS = 0.1 ! y-Startmittelpunkt vom Mittelpunkt beider
	! Wärmeübergangskoeffizienten in [mm]

	v1 = 0.5 ! Vorschubgeschwindigkeit der Düse im ersten
	! Streckenabschnitt in [mm s-1]

C	C v2 = 0.001E+3 ! Geschwindigkeit der Düse im zweiten
	! Streckenabschnitt in [m s-1]

C	C v3 = 0.0005E+3 ! Geschwindigkeit der Düse im dritten
	! Streckenabschnitt in [m s-1]

	x1 = 0.2 ! Streckenlänge des 1. Streckenabschnitts in [mm]
C	C x2 = 0.003E+3 ! Streckenlänge des 2. Streckenabschnitts in [mm]
C	C x3 = 0.004E+3 ! Streckenlänge des 3. Streckenabschnitts in [mm]
	time_Strecke1 = x1/v1 ! Zeit für die 1. Wegstrecke [s]
	timeperiod=time_Strecke1
C	C time_Strecke2 = x2/v2 ! Zeit für die 2. Wegstrecke [s]
C	C time_Strecke3 = x3/v3 ! Zeit die für die 3. Wegstrecke [s]
C	timeperiodstep1 = 1 ! Zeit des ersten Steps (in dem die
	! Anfangsbedingungen definiert werden) [s]
C	timeperiodstep2 = 1 ! Zeit des zweiten Steps (in dem nur der
	! Luftstrahl ohne Laser eingeschaltet wird
	! und sich nicht bewegt) in [s]
C	timeperiodstep3 = time_Strecke1 ! +time_Strecke2 + time_Strecke3
	! = Gesamtzeit des 3. Steps in dem der Laser und der
	! Luftstrahl eingeschalten sind [s]
	! !!!! muss berechnet und in Abaqus als "time period" eingetragen werden!!!!!

C----wenn der Mittelpunkt (xS) des gaußförmigen Wärmeübergangskoeff. die----
C--------Endposition (xEnde) erreicht hat, wird die Subroutine verlassen----
C-----------D.h. die Luftdüse wird "ausgeschaltet"------------------------
C	if(xS.gt.xEnde)then
C		RETURN
C		endif
		
C---------if-Schleife um den Luftstrommittelpunkt nach jedem Zeitinkrement--
C-----------------zu versetzen -> "Luftdüse bewegt sich"--------------------
C--„else if“-Bedingungen sind nur bei unterschiedl. Geschwindigkeiten nötig
	
C	if((TIME(2).gt.(timeperiodstep1 + timeperiodstep2)))then
C	if (TIME(1).le.time_Strecke1) then
	if (TIME(1).le.timeperiod) then
	xS = xS
	xS = v1*TIME(1) + xS
C	write(6,*) 'XS:',xS
	endif
C else if(TIME(1).lt.(time_Strecke1+time_Strecke2))then
C xS = v1*time_Strecke1 +xS
C xS = v2*(TIME(1)-time_Strecke1) + xS
C else if(TIME(1).lt.(time_Strecke1+time_Strecke2+time_Strecke3))then
C xS = v1*time_Strecke1 +v2*time_Strecke2
C xS = v3*(TIME(1)-(time_Strecke1+time_Strecke2)) + xS
C else
C xS = v1*time_Strecke1 +v2*time_Strecke2 +v3*time_Strecke3
C xS = v4*(TIME(1)-(time_Strecke1+time_Strecke2+time_Strecke3)) +xS
C	endif

	xS2 = xS + Abstand ! x-Startpunkt vom Mittelpunkt
	! des abgezogenen gaußförmigen Wärmeübergangsk.

C-------Berechnungen zum Wärmeübergangskoeffizienten------------------------
	xabstand = abs(COORDS(1)-xS) ! berechnet den x-Abstand vom
	! Knotenpunkt zum Mittelpunkt
	! des gaußförmigen W.k.

	xabstand2 = abs(COORDS(1)-xS2) ! berechnet den x-Abstand vom
	! Knotenpunkt zum Mittelpunkt des
	! abgezogenen gaußf. W.k.
	yabstand = abs(COORDS(2)-yS) ! berechnet den y-Abstand vom
	! Knotenpunkt zum Mittelpunkt
	! der Wärmeübergangskoeffizienten
	rabstand = sqrt(xabstand**2 + yabstand**2) ! berechnet den radialen
	! Abstand vom Knotenpunkt zum
	! Mittelpunkt des gaußf. W.k.
	rabstand2 = sqrt(xabstand2**2 + yabstand**2) ! berechnet den radialen
	! Abstand vom Knotenpunkt zum
	! Mittelpunkt des abgez. gaußf. W.k.

C------Hier wird der Wärmeübergangskoeffizient für jeden Knotenpunkt -------
C--------- je nach Position des Knotens berechnet---------------------------
C---if-Bedingung handelt die Streckung des Gaußprofils in eine Richtung ab--
C-- -> Temperaturfeld setzt sich aus Halbkreis und halber Ellipse zusammen--
	if((COORDS(1)-xS).le. 0)then
	gaussverteilung = z_Mittelpunktwaermekoeffizient*
     1 exp((-2)*xabstand**2/((Streckungsfaktor*Wo)**2)+
     2 (-2)*yabstand**2/Wo**2)-
     3 z_Mittelpunktwaermekoeffizient2*exp((-2)*rabstand2**2/Wo2**2)
       else
       gaussverteilung = z_Mittelpunktwaermekoeffizient*
     1 exp((-2)*rabstand**2/Wo**2)-
     2 z_Mittelpunktwaermekoeffizient2*exp((-2)*rabstand2**2/Wo2**2)
      endif

	if(KSTEP.gt.2)then !if(TIME(2).gt.1+time_Strecke1)then
	RETURN
	endif

      SINK=20 ! Fluidtemperatur [°C]
      H(1)=gaussverteilung ! Wärmeübergangskoeffizient [mW mm-2K-1]
      H(2)=0
      
C      	if(KINC.eq.1)then
C		write(6,*) H(1),';',xabstand,';',xS,';',KSTEP
C	endif

      RETURN
      END
C ###############################################################