#!/usr/bin/awk -f

function check_lang(prg, cfmt,x,y,lang,m1,m2) {
	cfmt = CONVFMT
	lang = ENVIRON["LANG"]
	CONVFMT = "%+.3e"
	x = 1.23
	y = "" x
	if("+1.230e+00"!=y) {
		# try to change language to "C" at run time
		ENVIRON["LANG"] = "C"
		x = 1.23
		y = "" x
		if("+1.230e+00"!=y) {
			m1 = ": numeric format of language '"
			m2 = "' not supported"
			print prg m1 lang m2 >> "/dev/stderr"
			m1 = ": set the environment variable LANG to C before running "
			print prg m1 prg >> "/dev/stderr"
			err = 1
			exit 1
		}
	}
	CONVFMT = cfmt
}

function print_usage(msg, m1) {
	if(msg!="") print "fam2dat: " msg >> "/dev/stderr"
	m1 = "[-h|-a] [OPT ...] SRC [OFILE1 [OFILE2]] [(-SP SPFILE) ...]"
	print "usage: fam2dat " m1 >> "/dev/stderr"
	print "OPT  : (-bl N) | (-p COL)" >> "/dev/stderr"
	print "SRC  : RUN | ( PFILE|-  QRUN|QFILE )" >> "/dev/stderr"
	print "SP is one of: LP BP HB PD TR UZ RO EP MX" >> "/dev/stderr"
	err = 1
	exit 1
}

function printb(msg,file) {
	if(append) {
		print msg >> file
	} else {
		if(lines[file]++ == 0) {
			print msg  > file
		} else {
			print msg >> file
		}
	}
	was_blank[file] = 0
}

function printa(msg,file) {
	if(msg == "") {
		if(!was_blank[file]) printb(msg,file)
		was_blank[file] = 1
	} else {
		printb(msg,file)
	}
}

# IBR  MTOT  ITP  LAB  NFPR  ISW  NTPL  NAR  NROWPR  NTST  NCOL  NPARX
#  1     2    3    4     5    6     7    8      9     10    11     12

function xtractsol(lab,prep, status,skip,i,M,N,LPP,lines,line,m1) {
	for(i in sol) delete sol[i]
	line_nos = 0
	
	for(status=(getline < sol_file); status==1; status=(getline < sol_file)) {
		if(NF>0 && $4==lab) break
		skip = $9
		for(i=0; i<skip; ++i) status=(getline < sol_file)
	}
	if( status != 1 ) {
		m1 = "error while looking for solution label "
		print "fam2dat: " m1 lab " in file " sol_file >> "/dev/stderr"
		exit 1
	}

	M=$7         # NTST*NCOL+1 number of mesh points
	N=$8         # NDIM+1      dimension of the problem
	LPP=int(N/7) #             lines per time point
	if(LPP*7<N) ++LPP
	lines=LPP
	line=""
	skip=$9

	for(status=(getline < sol_file); status==1; status=(getline < sol_file)) {
		if(NF>0) {
			--skip
			line=line $0
			if(--lines==0) {
				sol[line_nos++] = prep line
				line=""
				lines=LPP
				if(--M==0) break
			}
		}
	}
	if( status != 1 ) {
		m1 = "error while reading solution label "
		print "fam2dat: " m1 lab " from file " sol_file >> "/dev/stderr"
		exit 1
	}

	for(i=0; i<skip; ++i) status=(getline < sol_file)
	if( status != 1 ) {
		print "fam2dat: error while reading file " sol_file >> "/dev/stderr"
		exit 1
	}
}

function printsol(lab,prep,file, i) {
	if(lab != last_saved_lab) xtractsol(lab,prep)
	last_saved_lab = lab
	for(i=0; i<line_nos;    ++i) printb(sol[i], file)
	for(i=0; i<blank_lines; ++i) printb("",     file)
}

function exists(fname, status) {
	status = getline < fname
	close(fname)
	return ( status != -1 )
}

function run_name(run, pfname,qfname) {
	sol_file = ""
	if(run == "-") return run
	
	pfname = "data/p." run
	qfname = "data/q." run
	if(exists(pfname) && exists(qfname)) {
		sol_file = qfname
		return pfname
	}
	
	pfname = run
	if(exists(pfname)) return pfname
	
	print "fam2dat: run or file '" run "' does not exist" >> "/dev/stderr"
	err = 1
	exit 1
}

function sol_name(qfile, fname) {
	fname = "data/q." qfile
	if(exists(fname)) return fname
	
	fname = qfile
	if(exists(fname)) return fname
	
	print "fam2dat: qfile or file '" qfile "' does not exist" >> "/dev/stderr"
	err = 1
	exit 1
}

function next_sol( col) {
	if(curr_lab != 0) {
		last_lab  = curr_lab
		last_prep = curr_prep
		curr_lab  = 0
	} else {
		last_lab  = $4
		last_prep = ""
		for(col in cols) last_prep = last_prep $(col) " "
	}
}

function curr_sol( col) {
	curr_lab  = $4
	curr_prep = ""
	for(col in cols) curr_prep = curr_prep $(col) " "
}

BEGIN {
	check_lang("fam2dat")
	
	if(ARGC < 2 || ARGV[1] == "-h") print_usage()

	append      = 0
	blank_lines = 2
	curr_arg    = 1
	
	if(ARGC > curr_arg && ARGV[curr_arg]=="-a") {
		if(ARGC < 3) print_usage("too few arguments")
		append = 1
	}
	curr_arg += append

	for(;;) {
		if(ARGC > curr_arg && ARGV[curr_arg]=="-bl") {
			if(ARGC < 2+curr_arg+append) print_usage("too few arguments")
			blank_lines = ARGV[curr_arg+1]
			curr_arg += 2

		} else if(ARGC > curr_arg && ARGV[curr_arg]=="-p") {
			if(ARGC < 2+curr_arg+append) print_usage("too few arguments")
			cols[ARGV[curr_arg+1]+4]
			curr_arg += 2

		} else {
			break
		}
	}

	ARGV[1] = run_name(ARGV[curr_arg++])
	
	if(sol_file == "") {
		if(ARGC > curr_arg)
			sol_file = sol_name(ARGV[curr_arg++])
		else
			print_usage("too few arguments (QRUN missing)")
	}
	
	sb_out = "/dev/null"
	ub_out = "/dev/null"
	
	if( ARGC>curr_arg && substr(ARGV[curr_arg],1,1)!="-" )
		sb_out = ARGV[curr_arg++]

	ub_out = sb_out
	if( ARGC>curr_arg && substr(ARGV[curr_arg],1,1)!="-" )
		ub_out = ARGV[curr_arg++]
	
	show_stab = (sb_out != ub_out)
	curr_out  = sb_out

	lines[sb_out]     = 0
	lines[ub_out]     = 0
	was_blank[sb_out] = 1
	was_blank[ub_out] = 1
	
	sp_types[1]  = "-BP"
	sp_types[2]  = "-LP"
	sp_types[3]  = "-HB"
	sp_types[4]  = "-RO"
	sp_types[-4] = "-UZ"
	sp_types[5]  = "-LP"
	sp_types[6]  = "-BP"
	sp_types[7]  = "-PD"
	sp_types[8]  = "-TR"
	sp_types[9]  = "-EP"
	sp_types[-9] = "-MX"

	sp_args["-BP"]
	sp_args["-LP"]
	sp_args["-HB"]
	sp_args["-RO"]
	sp_args["-UZ"]
	sp_args["-PD"]
	sp_args["-TR"]
	sp_args["-EP"]
	sp_args["-MX"]

	for(arg=curr_arg; arg<ARGC; arg+=2) {
		if(ARGC>arg+1) {
			if(ARGV[arg] in sp_args) {
				sp_file = ARGV[arg+1]
				for(i in sp_types)
					if(sp_types[i]==ARGV[arg]) {
						sp_out[i]          = sp_file
						lines[sp_file]     = 0
						was_blank[sp_file] = 1
					}
			} else {
				print_usage("special point type " ARGV[arg] " not recognized")
			}
		} else {
			print_usage("syntax error at argument " ARGV[arg])
		}
	}
	
	bp_types[1] = "BP"
	bp_types[2] = "LP"
	bp_types[3] = "HB"
	bp_types[5] = "LP"
	bp_types[6] = "BP"
	bp_types[7] = "PD"
	bp_types[8] = "TR"

	state          =  0
	last_saved_lab = -1
	ARGC           =  2
}

# BR  PT  TY  LAB  PAR  L2-NORM  ...
#  1   2   3   4    5      6

              
$0 == "" { printa("",curr_out); state = 0; next }

$1 != 0 && $4 != 0 {
	pt_type = $3 % 10

	if(show_stab) {
		
		# begin stability splitting
		if       (state==0) {
			next_sol()
			curr_out = ($2<0) ? sb_out : ub_out
			state    = 1
			
		} else if(state==1) {
			stab = ($2<0) ? -1 : 1
			curr_out = (stab==1) ? ub_out : sb_out
			curr_sol()
			printsol(last_lab, last_prep, curr_out)
			next_sol()
			printsol(last_lab, last_prep, curr_out)
			state = 2
			
		} else if(stab*$2 > 0) {
			next_sol()
			printsol(last_lab, last_prep, curr_out)
			
		} else {
			pt_type = $3 % 10
			if(pt_type in bp_types) {
				curr_sol()
				printsol(curr_lab, curr_prep, curr_out)
				printa("", curr_out)
				stab = ($2<0) ? -1 : 1
				curr_out = (stab==1) ? ub_out : sb_out
				next_sol()
				printsol(last_lab, last_prep, curr_out)
			} else {
				printa("", curr_out)
				stab = ($2<0) ? -1 : 1
				curr_out = (stab==1) ? ub_out : sb_out
				curr_sol()
				printsol(last_lab, last_prep, curr_out)
				next_sol()
				printsol(last_lab, last_prep, curr_out)
 			}
			
		}
		#end stability splitting

	} else {
		next_sol()
		printsol(last_lab, last_prep, curr_out)
	}
	
	if(pt_type in sp_out) {
		printsol(last_lab, last_prep, sp_out[pt_type])
		printa("", sp_out[pt_type])
	}
}

END {
	if(state==1) printsol(last_lab, last_prep, curr_out)
	if(err == "") {
		for(file in lines) printa("",file)
	}
}
