def propagate_to_gain_goal(self, pulse_in, fiber, n_steps, power_goal = 1,
                              scalefactor_guess = None, powertol = 0.05):
        """
        Integrate over length of gain fiber such that the average output
        power is power_goal [W]. For this to work, fiber must have spectroscopic
        gain data from an amplifier model or measurement. If the approximate
        scalefactor needed to adjust the gain is known it can be passed as
        scalefactor_guess.\n This function returns a tuple of tuples:\n
        ((ys,AWs,ATs,pulse_out), scale_factor)
        """   
        if scalefactor_guess is not None:
                scalefactor = scalefactor_guess
        else:
                scalefactor = 1
        y, AW, AT, pulse_out = self.propagate(pulse_in,
                                                 fiber,
                                                 1,
                                                 output_power = scalefactor *\
                                                                 power_goal)         
        scalefactor_revised = (power_goal / (pulse_out.calc_epp()*pulse_out.frep))
        modeled_avg_power   = pulse_out.calc_epp()*pulse_out.frep                                                                 
        output_scale = scalefactor
        while abs(modeled_avg_power - power_goal) / power_goal > powertol:
            y, AW, AT, pulse_out = self.propagate(pulse_in, fiber,
                                                     1, output_power = 
                                                     power_goal * scalefactor_revised)
            modeled_avg_power_prev = modeled_avg_power
            modeled_avg_power = pulse_out.calc_epp()*pulse_out.frep
            slope = (modeled_avg_power - modeled_avg_power_prev) /\
                            (scalefactor_revised - scalefactor)                            
            yint = modeled_avg_power_prev - slope*scalefactor    
            # Before updating the scale factor, see if the new or old modeled
            # power is closer to the goal. When the loop ends, whichever is more
            # accurate is used for final iteration. This makes maximum use of 
            # each (computationally expensive) numeric integration
            if abs(modeled_avg_power-power_goal) < \
                        abs(modeled_avg_power_prev-power_goal):
                output_scale = scalefactor_revised
            else:
                output_scale = scalefactor
            # Update scalefactor and go to top of loop
            scalefactor = scalefactor_revised
            scalefactor_revised = (power_goal - yint)/slope                            
        print ('Updated final power:',modeled_avg_power,'W'         )
        return (self.propagate(pulse_in, fiber,
                              n_steps,
                              output_power = power_goal * output_scale),
                output_scale)