# wipe

puts "enter number of modes"
set numModes [gets stdin] 

# some settings for the response spectrum analysis
set tsTag 1; # use the timeSeries 1 as response spectrum function
puts "enter direction of seismic excitation 1=>X, 2=>Y, 3=>Z"
puts " "
set direction [gets stdin]; 

# define a 3D model

# the response spectrum function
set timeSeries_list_of_times_1 {0.0 0.06 0.1 0.12 0.18 0.24 0.3 0.36 0.4 0.42 \
								0.48 0.54 0.6 0.66 0.72 0.78 0.84 0.9 0.96 1.02 \
								1.08 1.14 1.2 1.26 1.32 1.38 1.44 1.5 1.56 1.62 \
								1.68 1.74 1.8 1.86 1.92 1.98 2.04 2.1 2.16 2.22 \
								2.28 2.34 2.4 2.46 2.52 2.58 2.64 2.7 2.76 2.82 \
								2.88 2.94 3.0 3.06 3.12 3.18 3.24 3.3 3.36 3.42 \
								3.48 3.54 3.6 3.66 3.72 3.78 3.84 3.9 3.96 4.02 \
								4.08 4.14 4.2 4.26 4.32 4.38 4.44 4.5 4.56 4.62 \
								4.68 4.74 4.8 4.86 4.92 4.98 5.04 5.1 5.16 5.22 \
								5.28 5.34 5.4 5.46 5.52 5.58 5.64 5.7 5.76 5.82 \
								5.88 5.94 6.0}
set timeSeries_list_of_values_1 {0.2 0.38 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.4762 \
								0.4167 0.3704 0.3333 0.303 0.2778 0.2564 0.2381 0.2222 0.2083 0.1961 \
								0.1852 0.1754 0.1667 0.1587 0.1515 0.1449 0.1389 0.1333 0.1282 0.1235 \
								0.119 0.1149 0.1111 0.1075 0.1042 0.101 0.098 0.0952 0.0926 0.0901 \
								0.0877 0.0855 0.0833 0.0813 0.0794 0.0775 0.0758 0.0741 0.0725 0.0709 \
								0.0694 0.068 0.0667 0.0641 0.0616 0.0593 0.0572 0.0551 0.0531 0.0513 \
								0.0495 0.0479 0.0463 0.0448 0.0434 0.042 0.0407 0.0394 0.0383 0.0371 \
								0.036 0.035 0.034 0.0331 0.0322 0.0313 0.0304 0.0296 0.0289 0.0281 \
								0.0274 0.0267 0.026 0.0254 0.0248 0.0242 0.0236 0.0231 0.0225 0.022 \
								0.0215 0.021 0.0206 0.0201 0.0197 0.0193 0.0189 0.0185 0.0181 0.0177 \
								0.0174 0.017 0.0167}
timeSeries Path 1 -time $timeSeries_list_of_times_1 -values $timeSeries_list_of_values_1 -factor 9.806

# define some analysis settings
constraints Transformation
numberer RCM
system UmfPack
test NormUnbalance 0.0001 10
algorithm Linear
integrator LoadControl 0.0
analysis Static

set eigs [eigen -genBandArpack $numModes]

# compute the modal properties
modalProperties -print -file "ModalReport.txt" -unorm

# define a recorder for the (use a higher precision otherwise the results
# won't match with those obtained from eleResponse)
set filename "eleforcesMode$direction.txt"
recorder Element -file $filename -closeOnWrite -precision 16 -eleRange 0 $m localForce

# currently we use same damping for each mode
set dmp [lrepeat [llength $eigs] 0.05]
# we don't want to scale some modes...
set scalf [lrepeat [llength $eigs] 1.0]
# CQC function
proc CQC {mu lambdas dmp scalf} {
	set u 0.0
	set ne [llength $lambdas]
	for {set i 0} {$i < $ne} {incr i} {
		for {set j 0} {$j < $ne} {incr j} {
			set di [lindex $dmp $i]
			set dj [lindex $dmp $j]
			set bij [expr [lindex $lambdas $i]/[lindex $lambdas $j]]
			set rho [expr \
				((8.0*sqrt($di*$dj)*($di+$bij*$dj)*($bij**(3.0/2.0))) / \
				(pow(1.0-$bij**2.0,2.0) + 4.0*$di*$dj*$bij*(1.0+pow($bij,2.0)) + \
				4.0*(pow($di,2.0) + pow($dj,2.0))*pow($bij,2.0)))]
			set u [expr $u + [lindex $scalf $i]*[lindex $mu $i] * [lindex $scalf $j]*[lindex $mu $j] * $rho]
		}
	}
	# it may happen, due to round-off erros for modes that are not excited in the requested direction,
	# that the final value of u is close to zero, but sligthly negative (for example -1e-20).
	# this is not good for sqrt!!!
	set u [expr max($u, 0.0)]
	return [expr sqrt($u)]
}

# ========================================================================
# TEST 01
# run a response spectrum analysis for each mode.
# then do modal combination in post-processing.
# ========================================================================
responseSpectrum $tsTag $direction

# read the input file for modal results.
# It should have N rows (N = number of modes) and M columns (M = 12 DOFs X Q elements).
# the output list will be 1 row x M columns. (1 = CQC combined result)
# finally we will create a file with Q x 12 (Q = number of elements, 12 = number of DOFs)

# real all input lines (1 line for each mode)
set f [open $filename "r"]
set lines [split [read $f] "\n"]
close $f
# save the total number of component just for the first line
set num_components 0
# convert it to a list of lists : M lists (1 for each component) of N entries (1 for each mode)
for {set i 0} {$i < [llength $lines]} {incr i} {
	# each line is a string with whitespace as separator
	set line [lindex $lines $i]
	if {[llength $line] > 0} {
		# convert the string into a list of strings
		set line [split $line " "]
		set num_components [llength $line]
		lset lines $i $line
	}
}
# now we can combine them into a single list
set combined_results {}
set modal_result [lrepeat $numModes 0.0]
# for each component....
for {set j 0} {$j < $num_components} {incr j} {
	# for each mode
	for {set i 0} {$i < $numModes} {incr i} {
		set cij [lindex [lindex $lines $i] $j]
		lset modal_result $i $cij
	}
	# combine this jth component
	set cqc_j [CQC $modal_result $eigs $dmp $scalf]
	# append it
	lappend combined_results $cqc_j
}

# now we can create a file with a line for each element
set cqc_filename "CQC_$filename"
set f [open $cqc_filename "w+"]
set counter 0
foreach item $combined_results {
	if {$counter > 0} {puts -nonewline $f " "}
	puts -nonewline $f $item
	set counter [expr $counter+1]
	if {$counter == 12} {
		puts -nonewline $f "\n"
		set counter 0
	}
}
close $f

# done
wipe
