# Frank Pulse version 1.0 
# by Christoffer Laustsen, cl@clin.au.dk, 2014

# Frank Pulse
proc Frankpulse {taup b1max phase N sweep} {
  global par

  set dtt [expr double($taup)/double($N)]
  set m [expr sqrt($N)]

  reset
  for {set i [expr -int($m)/2]} {$i < [expr int($m)/2]} {incr i} {
  	set k 0
  	for {set j [expr -int($m)/2]} {$j < [expr int($m)/2]} {incr j} {
		set t [expr $dtt*($k+0.5)]
		
		set amp [expr $b1max*(1.0-pow(cos(3.1415927*$t/($taup/$m)),$m))]
		set phi1 [expr 360/$m*fmod($sweep*$i*$j,$m)]
		
		pulse $dtt $amp $phi1
		set k [expr $k + 1]
  	}
  }
}

# Spin system spin = 1/2
spinsys {
    channels 1H
    nuclei   1H
    shift 1  0 0 0 0 0 0
}

# Parameters
par {
   spin_rate        0
   start_operator   I1z
   detect_operator  -I1z
   gamma_angles     1
   crystal_file     alpha0beta0
   sw               200000
   variable dt		2.0
   variable cp	   	1024
   variable sweep	0.512
   variable tpwr	2000
   variable N		256
   variable Tp		1024
   variable BW		4000
   np             	1
}

proc pulseq {} {
   global par
	
	reset
	Frankpulse $par(Tp) $par(tpwr) 0 $par(N) $par(sweep)
	
	acq
}

proc get_lims {SWH cp} {
  if { $cp <= 1} {
    set Res 0
  } else {
    set step [expr double($SWH)/($cp-1)]
    set Res {}
    for {set i 0} {$i < $cp} {incr i} {
      set shft [expr double($SWH)/2-$i*$step]
      lappend Res $shft
    }
  }
  return $Res
}

proc main {} {
   global par
  	
  set lims [get_lims $par(sw) $par(cp)]
  puts "Plot 1D Excitation effeincy: with SW: $par(sw) and points $par(cp)"
  
  set i 0
  set par(detect_operator) I1x
  puts "Start: $par(start_operator) Detect: $par(detect_operator)"
   set f [open $par(name)\_$par(sweep)\_$par(detect_operator) w]
        	foreach shft $lims {
	   incr i
	   set g [fsimpson [list [list shift_1_iso $shft]]]
           fscale $g -scale 2
            set ret [findex $g 1]
            puts $f "$shft $ret"
            puts "$shft $ret"
	    funload $g
        }
        close $f
  
    set i 0
  set par(detect_operator) I1y
  puts "Start: $par(start_operator) Detect: $par(detect_operator)"
   set f [open $par(name)\_$par(sweep)\_$par(detect_operator) w]
        	foreach shft $lims {
	   incr i
	   set g [fsimpson [list [list shift_1_iso $shft]]]
           fscale $g -scale 2
            set ret [findex $g 1]
            puts $f "$shft $ret"
            puts "$shft $ret"
	    funload $g
        }
        close $f
  
  set i 0
     set par(detect_operator) I1z
    puts "Start: $par(start_operator) Detect: $par(detect_operator)"
   	set f [open $par(name)\_$par(sweep)\_$par(detect_operator) w]
        	foreach shft $lims {
	   incr i
	   set g [fsimpson [list [list shift_1_iso $shft]]]
           fscale $g -scale 2
            set ret [findex $g 1]
            puts $f "$shft $ret"
            puts "$shft $ret"
	    funload $g
        }
        close $f
 }