#It is presumed that QuickField is already started, the problem coil.pbm is opened.


package require twapi;
twapi::import_commands;


#===============================================
#calculates force in QuickField
#===============================================
proc getQfForceX {x_position} {
#0. open model
#1. check plunger position, move plunger
#2. calculate force
	

# 0. open model
	set qf [comobj QuickField.Application];	#hook QuickField
    set pbm [[$qf Problems] Item 1];	#get the first opened Problem
    $pbm LoadModel; 					#load Model
    set mdl [$pbm Model];				#hook the Model
	
#1. check plunger x-position, the cooridate of vertex "R" on plunger
	set theVertex [[[[$mdl Shapes] Vertices] LabeledAs "R" "" ""] Item 1];
    set theXposition [[$theVertex Point] X];	#calculate x-coordinate of vertex "position"
    set theVector [$qf PointXY [expr $x_position - $theXposition] 0];	#calculate displacement vector
#plunger consists of 2 blocks: "insulation" (external), "plunger_core" (internal)
#if we move the external block then the internal block will be moved too	
    set theBlock [[[[$mdl Shapes] Blocks] LabeledAs "" "" "insulation"] Item 1];	#select "insulation" block	
    $theBlock Move 0 $theVector;		#move plunger. 0=qfShift, displacement



#2. Solve problem
    [$mdl Shapes] BuildMesh true false;	#build mesh
    $mdl Save;				#save model
    if {[$pbm Solved]!=1} {	#check if problem is solved
        $pbm SolveProblem;	#solve problem
		};
    $pbm AnalyzeResults;	#open results
    set res [$pbm Result];	#link to Result object
    

    set theFieldWindow [$res GetFieldWindow 1];	#hook Field window
    set theContour [$theFieldWindow Contour];	#build integration contour
#plunger consists of 2 blocks: "insulation" (external), "plunger_core" (internal)
    $theContour AddBlock1 "insulation";		#add to contour block "insulation"
    $theContour AddBlock1 "plunger_core";	#add to contour block "plunger_core"

#calculate magnetic force integral
    set theForce [[$res GetIntegral 15 $theContour] Value]; #15=qfInt_MaxwellForce
	return [$theForce X]; 	#get x-component of force
	}

	
	
	

#===============================================
# main programm starts here
#===============================================

#relay physical parameters
set x_out 0.01;       		# pull out position. m
set x_in 0.006;       		# pull in position, m
set m_mass 0.0045;     		# plunger weight, kg
set k 4.0;            		# spring constant, N/m
set x_spring_free 0.015;  	# spring free would be position, m


#initial conditions
set t_time 0.0;    			# time, s
set x_position $x_out;		# plunger position, m
set v_speed 0.0;			# plunger speed, m/s


set dt_time_step 0.005;		#integration time step, s
set max_time 0.2;			# maximal integration time


puts "Time,s	Position,m	Speed,m/s";

while { $t_time < $max_time && $x_position > $x_in } {
	set theForce [getQfForceX $x_position];			# electromagnetic force
    set theForce [expr $theForce + ($x_spring_free - $x_position) * $k];	# add spring force
    set v_speed [expr $v_speed +  ($theForce / $m_mass) * $dt_time_step]; 	# calculate acceleration and speed
    set x_position [expr $x_position + $v_speed * $dt_time_step];			# calculate plunger position
	set t_time [expr $t_time + $dt_time_step];		#time increment
	puts "$t_time	$x_position	$v_speed";			#output results
	};